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Abstract 


This  dissertation  develops  a  new  wavelet  design  technique  that  produces  a  wavelet  that  matches  a  desired 
signal  in  the  least  squares  sense.  The  Wavelet  Transform  has  become  very  popular  in  signal  and  image 
processing  over  the  last  6  years  because  it  is  a  linear  transform  with  an  infinite  number  of  possible  basis 
junctions  that  provides  localization  in  both  time  (space)  and  frequency  (spatial  frequency). 

The  Wavelet  Transform  is  very  similar  to  the  matched  filter  problem,  where  the  wavelet  acts  as  a  zero 
mean  matched  filter.  In  pattern  recognition  applications  where  the  output  of  the  Wavelet  Transform  is 
to  be  maximized,  it  is  necessary  to  use  wavelets  that  are  specifically  matched  to  the  signal  of  interest. 
Most  current  wavelet  design  techniques,  however,  do  not  design  the  wavelet  directly,  but  rather,  build  a 
composite  wavelet  from  a  library  of  previously  designed  wavelets,  modify  the  bases  in  an  existing  mul¬ 
tiresolution  analysis  or  design  a  multiresolution  analysis  that  is  generated  by  a  scaling  Junction  which 
has  a  specific  corresponding  wavelet.  In  this  dissertation,  an  algorithm  for  finding  both  symmetric  and 
asymmetric  matched  wavelets  is  developed.  It  will  be  shown  that  under  certain  conditions,  the  matched 
wavelets  generate  an  orthonormal  basis  of  the  Hilbert  space  containing  all  finite  energy  signals.  The 
matched  orthonormal  wavelets  give  rise  to  a  pair  of  Quadrature  Mirror  Filters  (QMF)  that  can  be  used 
in  the  fast  Discrete  Wavelet  Transform.  It  will  also  be  shown  that  as  the  conditions  are  relaxed,  the 
algorithm  produces  dyadic  wavelets  which  when  used  in  the  Wavelet  Transform  provides  significant  re¬ 
dundancy  in  the  transform  domain. 

Finally,  this  dissertation  develops  a  shift,  scale  and  rotation  invariant  technique  for  detecting  an 
object  in  an  image  using  the  Wavelet  Radon  Transform  (WRT)  and  matched  wavelets.  The  detection 
algorithm  consists  of  two  levels.  The  first  level  detects  the  location,  rotation  and  scale  of  the  object, 
while  the  second  level  detects  the  fine  details  in  the  object.  Each  step  of  the  wavelet  matching  algorithm 
and  the  object  detection  algorithm  is  demonstrated  with  specific  examples. 
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Chapter  1 


Introduction 


For  many  years,  the  primary  linear  mathematical  analysis  tool  used  to  transform  signal  information  was 
the  Fourier  Transform.  However,  the  end  of  the  1980’s  saw  the  development  of  an  alternative  mathe¬ 
matical  framework,  called  the  wavelet  transform,  with  applications  in  signal  and  image  analysis  [30]. 
While  much  of  the  theory  associated  with  it  is  not  new  and  can  be  described  in  a  Hilbert  space  setting, 
it  did  provide  a  consolidated  framework  for  a  number  of  previously  diverse  disciplines,  like  multireso¬ 
lution  analysis  used  in  computer  vision,  subband  coding  developed  for  speech  and  image  compression, 
and  orthonormal  basis  expansions  developed  in  applied  mathematics  [30]. 

The  Fourier  Transform  of  a  signal,  f{x),  given  by 

/OO 

f(x)e~^‘^^dx,  (1.1) 

-OO 

is  a  projection  operation  onto  the  basis  formed  by  dilating  the  complex  exponential,  The  Wavelet 
Transform,  given  by 

W4a,b)  =  (^)  dx,  (1.2) 

is  a  projection  operation  onto  the  basis  formed  by  dilating  and  shifting  a  “mother”  wavelet,  '>p{x),  which 
is  zero  mean  and  decays  very  rapidly  in  x.  These  properties  of  the  wavelet  provide  an  advantage  over 
Fourier  analysis  because  the  Wavelet  Transform  can  achieve  localization  in  both  time  and  frequency  or 
in  the  case  of  images,  space  and  spatial  frequency  [10].  That  is  to  say  that  the  Wavelet  Transform  can 
provide  information  about  the  frequency  content  of  a  particular  group  of  pixels.  Another  advantage  of 
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wavelet  analysis  over  Fourier  analysis  is  the  flexibility  in  the  transform  operator  used  for  decomposition 
or  signal  representation.  Fourier  methods  are  locked  into  bases  that  are  based  on  the  complex  exponential 
kernel,  e*®.  Wavelet  transforms  can  have  many  different  operators.  They  may  or  may  not  be  orthogonal, 
and  they  may  or  may  not  form  bases.  Even  in  the  case  of  orthonormal  bases,  there  are  still  infinitely  many 
wavelets  that  satisfy  the  appropriate  conditions.  Because  of  this  flexibility,  it  is  possible,  and  wise,  to 
choose  or  design  a  wavelet  that  is  matched  to  the  application. 

Many  wavelet  design  techiques  have  been  developed  since  the  early  theoretical  work  of  the  wavelet 
pioneers.  However,  most  of  the  techniques  do  not  directly  match  a  wavelet  to  a  signal  of  interest.  Many 
build  adaptive  wavelets  from  existing  wavelets  and  their  effectiveness  is  dependent  on  the  effectiveness 
of  the  existing  wavelets  used.  Others  impose  requirements  on  the  wavelets  to  ensure  some  degree  of 
smoothness  or  differentiability,  which  guarantees  the  wavelet  to  exhibit  certain  properties,  but  does  not 
guarantee  anything  with  regards  to  its  shape. 

Using  adaptive  wavelets  for  image  pattern  recognition  is  very  attractive  because  of  the  scale  (or 
zoom)  insensitivity  of  the  wavelet  transform.  For  instance,  if  a  wavelet  can  be  constructed  that  matches 
a  pattern  of  interest  in  an  image,  then  the  peak  of  the  wavelet  transform  (or  lack  thereof)  will  indicate 
whether  the  pattern  is  present  (or  not)  and  due  to  the  localization  properties  of  the  wavelet,  where  the 
pattern  is  located.  Applying  the  Continous  Wavelet  Transform  (CWT)  to  the  image  would  be  too  data 
intensive  because  of  the  redundancy  contained  in  the  CWT.  The  2-D  Discrete  Wavelet  transform  is  much 
faster  because  it  processes  the  projection  coefficients  only.  However,  it  is  not  shift  or  rotation  invariant, 
making  its  application  to  pattern  recognition  of  limited  use. 

These  two  problems,  wavelet  matching  and  application  to  pattern  recognition  are  addressed  in  this 
dissertation.  First,  a  wavelet  design  algorithm  is  developed  that  takes  any  1-D  signal  as  input  and  finds 
the  wavelet  that  comes  closest  to  it  in  the  least  squares  sense.  The  algorithm  assumes  that  the  wavelet 
is  bandlimited.  The  bandlimit  conditions  can  be  set  such  that  the  resultant  wavelet  is  an  orthonormal 
basis  of  the  function  space  or  they  can  be  relaxed  or  widened,  so  that  the  algorithm  generates 

a  “dyadic”  wavelet.  Second,  a  target  detection  and  identification  algorithm  is  developed  based  on  the 
Radon  Transform,  and  the  Continuous  Wavelet  Transform  using  matched  wavelets.  The  target  detec- 
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tion  algorithm  is  shift,  scale  and  rotation  invariant.  The  Radon  Transform  is  used  to  reduce  the  data 
burden  on  the  CWT  and  to  provide  some  degree  of  rotation  invariance.  Using  matched  wavelets  in  the 
CWT  provides  scale  and  shift  invariance  as  well  as  optimal  amplitude  detection.  The  matched  wavelet 
algorithm  in  conjunction  with  the  target  detection  algorithm  provide  a  thorough  framework  from  which 
more  detailed  target  classification  algorithms  could  be  developed. 

This  dissertation  provides  an  in-depth  review  of  wavelet  theory,  first  in  Chapter  2  with  a  short  review 
of  Hilbert  spaces  and  a  comparison  between  Fourier  and  Wavelet  analysis.  Chapter  3  provides  an  histor¬ 
ical  chronology  of  the  development  of  wavelet  theory  from  the  Continuous  Wavelet  Transform  of  Morlet 
and  Grossman  to  multiresolution  analyses  and  the  Discrete  Wavelet  Transform  of  Mallat.  Several  current 
wavelet  design  techniques  are  presented  in  Chapter  4  as  motivation  for  the  matching  algorithm  devel¬ 
oped  in  Chapter  5.  Chapter  6  gives  several  examples  of  the  wavelet  matching  algorithm  to  demonstrate 
its  performance  and  to  assist  anyone  tiying  to  implement  the  algorithm  themselves.  After  a  summary 
of  the  Radon  Transform  and  image  reconstruction  using  backprojection.  Chapter  7  provides  the  step  by 
step  development  of  the  target  detection  and  identification  algorithm  taking  advantage  of  the  matched 
wavelet  algorithm  of  Chapter  5.  Chapter  8  provides  a  summary  of  the  research  results  contained  in  this 
dissertation,  and  the  Appendix  contains  proofs  to  many  of  the  theorems. 
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Chapter  2 


Background 


This  chapter  contains  brief  background  material  that  will  be  helpful  in  understanding  the  development 
of  wavelet  analysis.  The  first  section  contains  definitions  of  terms  used  throughout  the  dissertation,  and 
an  introduction  to  Hilbert  spaces,  the  projection  theorem  and  its  application  to  signal  representations 
in  Hilbert  spaces.  The  Fourier  series  is  shown  to  be  a  direct  application  of  the  projection  theorem  on 
the  Hilbert  space  of  all  periodic,  finite  energy  signals.  The  limitations  of  Fourier  analysis  on  real  world 
signals  are  indentified  followed  by  the  introduction  of  the  Short  Time  Fourier  Transform  (STFT).  The 
last  section  introduces  wavelet  analysis  as  a  recent  solution  to  the  limitations  of  Fourier  analysis  and  the 
STFT.  More  detailed  development  of  wavelet  theory  will  be  provided  in  Chapter  3. 


2.1  Hilbert  Spaces 


The  Hilbert  space  is  a  complete,  linear  vector  space  with  a  norm  ||  •  ||  and  an  inner  product  (•,  •)  defined 
[22].  There  are  two  Hilbert  spaces  used  in  wavelet  analysis,  L'^{U),  and,  where  3?  is  the  set  of 

all  real  numbers  and  Z  is  the  set  of  all  integers. 


Definition  1  The  Hilbert  space  consists  of  complex-valued,  measurable  functions,  x(t),  on  the 
real  line,  3?,  where 


/oo 

\x{t)\‘^dt  <  oo 

-oo 


(2.1) 
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and  the  integral  is  the  Lebesque  integral.  The  norm  of  x  e  L'^{W)  is  defined  as 


ikii = (y  \x{t)\^d^^ . 

The  inner  product  ofx,  y  €  1-^(3?)  is  defined  as 

/OO  _ 

x{t)y{t)dt 

-00 


where  y{t)  is  the  complex  conjugate  ofy{t). 


(2.2) 


(2.3) 


Definition  2  The  Hilbert  space  consists  of  all  complex  valued  sequences  of  scalars 

X  =  {...  ...}  for  which 

OO 

Y.  1^*1^  <  OO  (2.4) 

i=— OO 

for  all  i  ^  Zy  the  integer  number  line.  The  norm  ofx  is  defined  as 

•  (2.5) 

The  inner  product  ofx  =  {. . .  5?7_i,  •  •  •}  z/q,  . . .}  is  defined  as 


OO 

Y  (2.6) 

i=—oo 

A  set  of  vectors  in  a  Hilbert  space,  Xi  €  H,  is  said  to  be  orthonormal  if  Xi  ±  xj  for  i  7^  j  for  all  Xi  e  H 
and  if  each  vector  has  unit  norm,  that  is. 


0 


for  i  =  j 
for  i  /  j 


(2.7) 


The  projection  theorem  is  a  fundamental  theorem  used  in  vector  space  analysis  and  is  used  to  formulate 
the  mathematics  of  both  Fourier  and  wavelet  decompositions  for  functions  in  Hilbert  spaces. 


Theorem  1  Projection  Theorem  Let  H  be  a  Hilbert  space  and  M  a  closed  subspace  of  H.  Corre¬ 
sponding  to  any  vector  x  e  H,  there  is  a  unique  vector  mo  e  M  such  that  ||a;  -  mo||  <  ||x  -  m||/or 
all  m  €  M.  Furthermore,  a  necessary  and  sufficient  condition  that  mo  €  M  be  the  unique  minimizing 
vector  is  that  x  -  mo  be  orthogonal  to  M  [22  J. 
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The  Projection  Theorem  states  that  if  we  want  to  construct  a  vector  mo  €  M,  a  subspace  of  H,  then  we 
can  do  so  uniquely.  Furthermore,  the  vector  that  best  approximates  x  will  be  the  one  that  minimizes  the 
norm  of  the  error  vector,  |(a;  —  mo  ||,  and  the  error  vector,  x  —  ttiq,  will  be  orthogonal  to  the  approximation 
subspace,  M  (Figure  2.1)  [22]. 

Suppose  we  want  to  construct  an  approximation  oixeH  from  a  set  of  vectors  y  =  {yi ,  y2 ,  •  •  • ,  l/n  }  e 
M ,  where  M  is  a  closed  subspace  of  H  and  y  is  a  basis  of  M .  Then,  the  approximation,  or  projection 


Figure  2. 1 :  Orthonormal  projection 
of  X  onto  M  is  a  linear  combination  of  y^ 

n 

{Pmx)  =  am  (2.8) 

i=l 

where  Pm  is  the  projection  operator  onto  the  subspace  M  [22].  From  the  Projection  Theorem,  the  best 
approximation  of  x  will  be  the  one  that  minimizes  \\x  -  Pmx\\  and  x  -  Pmx  will  be  orthogonal  to  all 
elements  of  y,  that  is, 

{x  -  PMX,yi)  =  0.  (2.9) 
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Substituting  (2.8)  into  (2.9)  gives 

n 

{x-Y^ajVj.yi)  =  Q.  (2.10) 

Since  the  Hilbert  space  is  linear  and  the  inner  product  is  a  linear  operator,  (2.10)  can  be  rewritten  as 

n 

{^>yi)  =  '^oij{yj,yi).  (2.11) 

j=l 

If  y  €  M  is  an  orthonormal  basis  of  M,  then  {yi,yj)  =  0  for  i  ^  j  and  ||yi||  =  1.  Equation  (2.11) 
reduces  to 

{x,yi)  =  ai  (2.12) 

making  (2.8) 

n 

{PMx)  =  Y^{x,yi)yi.  (2.13) 

i=l 

Equation  (2.13)  gives  the  expression  for  the  projection  of  the  vector  x  onto  an  orthonormal  basis  y  of 
M  [22]. 

2.2  Fourier  Analysis 

Signal  and  image  analysis  has  long  benefited  from  Fourier  analysis,  where  the  set  of  functions, 
forms  an  orthonormal  basis  of  L^(0,T),  the  Hilbert  space  consisting  of  all  square  integrable  functions 
defined  on  the  interval  [0,T].  All  periodic  functions,  with  period  T,  also  reside  in  L^(0,r),  since  they 
can  be  completely  described  by  one  period  on  the  interval  [0,T].  Since  L'^{0,T)  is  a  Hilbert  space,  it 
has  an  inner  product  and  a  norm  [10]  defined  as 

=  f{x)g{x)dx  (2.14) 

\\f\\  =  {fJ)=(^^l^\f{x)\^dx'^  \  (2.15) 

Since  is  an  orthonormal  basis  of  f{x)  e  L‘^{0,T)  can  be  represented  using  (2.13) 

and  (2.12), 

OO 

/(x)=  ^  ckc^^  (2.16) 

fc  =  ~00 
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(2.17) 


c*  =  f{x)e  dx. 

The  function,  /,  is  decomposed  into  the  sum  of  infinitely  many  mutually  orthogonal  components, 
9kix)  =  where  orthogonality  means: 

{9n,gm)=0,  for  all  m  7^  n.  (2.18) 

That  (2.18)  holds  is  a  consequence  of  the  fact  that: 

Wk{x)  =  6*2-*=-/^  =  {...,  -1,0, 1, .. .}  (2.19) 

is  an  orthonormal  basis  of  7)2(0,  T)  [10].  It  is  important  to  note  further  that  the  orthonormal  basis,  Wk, 
is  formed  by  dilating  (known  as  integral  dilation)  a  single  function,  Chui  [10]  emphasizes  the 
significance  of  these  facts  in  the  introduction  to  his  book: 

Let  us  summarize  this  remarkable  fact  by  saying  that  every  2n-periodic  square  integrable 
function  is  generated  by  a  "superposition  ”  of  integral  dilations  of  the  basic  Junction  w{x)  = 

Chui  develops  his  section  on  Fourier  analysis  assuming  a  period  of  2;r,  thereby  giving  rise  to  his  mention 
of  27r -periodic  functions  and  the  basic  function  e*®. 

The  decomposition  in  (2.16)  and  (2.17)  is  known  as  the  Fourier  series  expansion  of  /(a;)[17].  The 
basis,  is  a  complex  exponential  with  fundamental  frequency,  =  l/T.  Integer  multiples  of 

the  fundamental  frequency  are  called  harmonics.  So,  a  periodic  function  of  period  T,  can  be  perfectly 
represented  by  the  sum  of  a  complex  exponential  at  a  fundamental  frequency  and  its  harmonics.  The 
weights  of  these  complex  exponential  components,  Cfc,  constitute  the  frequency  spectrum  off.  It  is  often 
called  the  line  spectrum  because  it  is  discrete. 

Square  integrable,  non-periodic,  functions  whose  domain  is  the  real  number  line,  5ft,  constitute  the 
Hilbert  space,  i)2(5ft),  with  its  inner  product  and  norm  defined  as: 


POO 

{f^9)=  f{x)g{x)dx 

*/ — OO 

(2.20) 

(2.21) 

8 


The  Fourier  series  of  / (x)  E  L^(St)  does  not  exist  since  /  is  non-periodic  and  its  domain  is  the  entire 
real  line,  3?.  The  frequency  spectrum  of  /  can  still  be  determined  by  taking  its  inner  product  with  the 
function,  where  ^  is  the  continuous  frequency  variabte[17].  Letting  cj  =  27t^  gives 

/OO 

f{x)e-^^^dx.  (2.22) 

-CX) 

F{ixi)  is  the  Fourier  Transform  of  f{x)  and  is  at  least  piecewise  continuous  [1 7],  The  Fourier  Transform 
can  be  interpreted  in  a  manner  similar  to  the  Fourier  series.  A  non-periodic,  square  integrable  function, 
f{x)  E  L^(3i)  can  be  represented  by  the  integral  sum  of  complex  exponentials  with  weights  given  by 
the  frequency  spectrum,  F(w), 

1  C°° 

fi^)  =  TT  j  F{u))e^'^du.  (2.23) 

ZTT  7—00 

The  Fourier  Series  and  Fourier  Transform  are  powerful  tools  for  determining  the  frequency  content 
of  discrete  and  continuous  signals.  However,  in  order  to  study  the  spectral  behavior  of  analog  signals 
from  its  Fourier  Transform,  full  knowledge  of  the  signal  in  the  time  or  space  domain  must  be  obtained[  1 0]. 
In  addition,  if  a  signal  is  altered  in  a  small  neighborhood  of  some  time  instant,  then  the  entire  spectrum 
is  affected.  If  a  signal  of  some  frequency  exists  for  a  finite  period  of  time,  the  Fourier  Transform  does 
not  give  visibility  into  where  in  time  the  signal  occurred.  Being  able  to  determine  where  in  time  a  sig¬ 
nal  of  some  frequency  occurs  is  called  time  localization.  Being  able  to  determine  the  frequency  content 
of  a  signal  at  a  particular  frequency  is  called  frequency  localization.  The  Fourier  Transform  provides 
excellent  frequency  localization,  but  poor  time  localization  [10]. 

In  order  to  obtain  localization  in  both  time  and  frequency,  an  additional  parameter  must  be  added. 
Gabor  [16]  was  the  first  to  adapt  Fourier  analysis  to  include  a  modulation  window,  g{x)  [30].  Gabor’s 
Short-Time  Fourier  Transform  (STFT)  [30],  5(r,  w),  is  given  as 


/OO  _ 

fix)g{x  -  T)e~^'^dx. 

-OO 


(2.24) 


The  signal  is  weighted  by  a  finite  duration  window,  g{x),  prior  to  taking  the  Fourier  Transform  (Figure 
2.2).  The  additional  parameter,  r,  is  the  translation  parameter  of  the  window  and  it  provides  the  addi¬ 
tional  dimension  needed  to  obtain  localization  in  both  time  and  frequency.  If  a  frequency  pulse  or  burst 
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Figure  2.2:  Short-Time  Fourier  Transform 
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is  present  in  the  signal  for  a  finite  duration,  it  will  show  up  in  at  values  of  r  corresponding  to 

the  location  of  the  frequency  pulse.  The  function,  S{t,  u),  provides  a  two  dimensional  map  of  time  and 
frequency  (Figure  2.3)  [30].  By  fixing  frequency  at  some  value,  wq,  S(t,u;o)  gives  a  1-D  function  indi- 


At 

◄ - ► 

Ao) 

Figure  2.3:  Time- frequency  map  -  Short  Time  Fourier  Transform 


eating  where  in  time  wq  existed.  Alternatively,  by  fixing  r  =  tq,  S(to,  u)  gives  the  Fourier  Transform 
for  that  time  slice  of  / .  The  resolution  in  both  time  and  frequency  depend  on  the  window  function  g{t). 
Time  and  frequency  resolutions  are  traded  off  according  to  the  uncertainty  principle  [10,  30] 


Time-bandwidth  product  =  AtAu  >  - 

2 


(2.25) 


where  At  and  Acu  are  the  rms  measures  of  time  width  and  bandwidth  given  by 


At 


{/f^  (a:  -  E{g))  \g{x)\^  dx^ 

M 


Au; 


_  {/f°^(a;-,E(G))|G(cc;)|2du;} 

IIGII 


(2.26) 


(2.27) 
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G(u;)  is  the  Fourier  Transform  of  g(x)  and  the  operator,  B(-)  is  given  as 


ll/IP 


(2.28) 


For  this  reason,  a  gaussian  window  is  often  used  since  it  satisfies  the  equality  of  (2.25)  [30].  However, 
modulating  the  signal  with  a  gaussian  window  correlates  the  Fourier  coefficients  and  therefore  destroys 
orthonormality.  Furthermore,  the  window  function,  g(x),  with  its  corresponding  frequency  spectrum, 
G(u;),  sets  both  the  time  and  frequency  resolutions  for  the  entire  time-frequency  plane.  This  is  a  disad¬ 
vantage  for  analyzing  signals  with  both  high  and  low  frequencies  [10].  In  order  to  properly  represent  a 
signal  of  some  fundamental  frequency,  the  window  must  contain  one  or  more  periods  of  that  signal.  Low 
frequency  signals,  therefore,  require  long  time  windows,  which  corresponds  to  high  resolution  in  the  fre¬ 
quency  domain.  High  flrequency  signals  require  a  small  time  window  in  order  to  capture  one  or  more 
periods.  The  small  time  window  corresponds  to  low  resolution  in  the  frequency  domain.  The  STFT’s 
window  widths  are  fixed  in  both  time  and  frequency,  as  illustrated  in  Figure  2.3,  and  therefore,  cannot 
effectively  analyze  signals  containing  both  high  and  low  frequencies  [10]. 

The  obvious  solution  to  Gabor’s  Short-Time  Fourier  Transform  is  an  orthonormal  basis  of  L^(5i) 
that  provides  both  time  and  frequency  localization  for  signals  with  high  and  low  frequencies.  Wavelet 
analysis  provides  just  such  a  solution. 


2.3  Wavelet  Analysis 


The  “basic  wavelet”  or  “mother  wavelet”  is  a  function,  i/)(x)  G  L^(3fJ),  whose  Fourier  transform,  ’J'(a;), 
satisfies  the  admissibility  condition  [10,  12,  14,  23,  27]: 


Since  ^(x)  G  L^(3t),  then 


|w| 


du)  <  00. 


dx  <  oo 


(2.29) 


(2.30) 
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must  be  trae,  which  means  that  ^(a;)  must  have,  effectively,  finite  support.  Assuming  ^(cj)  is  continu¬ 
ous,  (2.29)  and  (2.30)  imply 


T  f°° 

’i'(O)  =  0,  /  il’{x)dx  =  0. 

J~oo 


(2.31) 


This  is  the  reason  that  'tp  is  called  a  “wavelet”  [10].  In  order  for  ip  to  have  an  average  of  0  (2.3 1 ),  it  must 
be  wave-like  in  nature.  In  order  for  it  to  be  in  T2(3fi),  it  must  be  effectively  finite  or  of  short  duration  as 
in  (2.30),  hence  the  term  “small  wave”  or  “wavelet.”  Equation  (2.31)  also  indicates  that  ip(x)  has  the 
characteristics  of  a  bandpass  filter. 

If  we  let 


^a,b  =  |a|  2^P{- - 

a 

then  the  continuous  wavelet  transform  (CWT)  [10, 12,  30]  is  defined  as 


(2.32) 


1  roo  rp  _ _  A 

Wf{a,h)^{f,ipafi)  =  I  f{x)ip{ - )dx  (2.33) 

where  / (x)  G  1/^(3?),  a;  6  G  3?,  and  a  /  0.  The  parameters,  a  and  b,  are  the  scale  and  shift  parameters, 
respectively,  and  they  provide  localization  in  both  the  time  and  frequency  domains.  The  parameter  b 
centers  the  wavelet  at  f  =  6  and  a  scales  the  wavelet  function,  ip,  on  the  a:-axis.  The  CWT  is  analogous 
to  the  Fourier  Transform  where  frequency  has  been  replaced  by  scale,  since  a  gives  insight  into  the  fre¬ 
quency  content  of  a  passband,  not  a  single  frequency.  The  bandwidth  of  the  passband  changes  with  o 
[10,  30]  such  that 

Au 

- =  K.  (2.34) 

CJ 

This  characteristic  is  common  in  communication  theory  and  is  called  “constant  Q”  (Figure  2.4)  [10, 30]. 
Now,  the  windowing  function  is  effectively  the  support  of  the  wavelet,  ip,  which  changes  with  a.  It  is 
not  fixed  as  in  the  case  of  the  STFT,  and  can  therefore  provide  both  time  and  frequency  localization  of 
a  signal  containing  both  high  and  low  frequencies. 
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Chapter  3 


Wavelet  Theory  Development 


In  this  chapter,  the  history  of  the  development  of  wavelet  theory  will  be  presented  from  the  original  work 
of  Morlet  and  Grossman  through  the  algorithm  development  of  Stephane  Mallat.  The  development  of 
wavelet  theory  was  motivated  by  the  need  to  analyze  a  finite  energy  signal  with  a  single,  finite  energy 
function,  called  a  wavelet,  dilated  and  shifted  by  real  parameters.  It  was  shown  that  the  condition  on  the 
analyzing  function  was  simple  and  easily  met.  As  the  constraints  on  the  dilation  and  shift  parameters 
tightened,  for  instance,  confinement  to  2L,  the  conditions  on  the  analyzing  function  also  increased.  It 
was  shown  by  Meyer  that  if  the  dilation  parameter  values  were  constrained  to  be  powers  of  2  and  the 
shift  parameter  values  to  be  integer  multiples  of  powers  of  2,  then  a  wavelet  could  be  found  such  that 
its  shifts  and  dilates  formed  an  orthonormal  basis  of  L^(5R).  This  chapter  goes  through  in  chronological 
order  the  development  of  each  “class”  of  wavelet  and  the  conditions  associated  with  each,  culminating 
with  orthonormal  bases  and  multiresolution  analyses. 

3.1  Continuous  Wavelet  Transform 

J.  Morlet,  a  French  geophysicist,  first  proposed  the  use  of  “wavelets  of  constant  shape”  for  analyzing 
seismic  data.  His  reference  to  constant  shape  was  intended  to  contrast  these  new  functions  with  the 
Short  Time  Fourier  Transform  (STFT),  which  are  not  of  constant  shape  [14],  A.  Grossman,  a  French 
theoretical  physicist,  showed  that  Morlet’s  function,  dilated  and  shifted  by  real  parameters,  generated 
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a  square  integrable  representation  of  functions  in  L^QR.)  and  this  representation  was  referred  to  as  the 
wavelet  transform  [14,  30],  The  wavelet  transform  is  defined  as  follows: 


Definition  3  (Continuous  Wavelet  IVansform)  [lOJlf-ip  e  satisfies  the  “admissibility"  con 
dition: 


~  /  — i — i — ^ 

J  —  OO  1^1 

where  is  the  Fourier  transform  of  V',  then  xj}  is  called  a  “basic  wavelet", 
wavelet  the  continuous  wavelet  transform  (CWT)  on  is  defined  by 


(3.1) 


Relative  to  every  basic 


Wf{a,b)  =  {f,i)a,b) 


(3.2) 


for  f  6  L2(3?),  and  a,  6  €  3?  with  o  7^  0. 


The  wavelet  operator, 

'il’a,b{x)  =  \a\~^'ip  >  (3.3) 

is  a  dilated  and  shifted  version  of  a  single  function,  'ijj{x),  sometimes  referred  to  as  the  mother  wavelet, 
and  it  maps  a  one  dimensional  signal  into  a  two  dimensional  analysis  domain,  that  is,  scale,  a,  and  shift, 
b.  This  mapping  provides  visibility  into  the  frequency  content  of  a  signal  through  the  scale  parameter, 
a,  and  time  localization  through  the  shift  parameter,  6,  a  significant  advantage  over  standard  Fourier 
analysis.  Figure  3. 1  shows  the  effect  of  the  dilation  and  shift  parameters,  a  and  b  on  the  wavelet  operator. 
The  wavelet  shown,  used  by  Morlet  andGrossman,  is  a  cosine  function  modulated  by  a  gaussian  window. 
Notice  that  as  a  increases,  the  wavelet  gets  wider  and  shorter  in  height  and  when  a  decreases,  the  wavelet 
gets  thinner  and  taller.  This  effect  provides  a  zoom-in,  zoom-out  capability  in  the  frequency  domain. 
The  bandwidth  of  the  wavelet  for  small  a  is  large  and  the  bandwidth  for  large  a  is  small.  Of  course, 
changes  in  b  shift  the  wavelet  up  and  down  the  a;-axis  giving  full  two  dimensional  latitude  in  both  time 
and  frequency.  The  admissibility  condition  (3.1)  on  implies 

'F 

’^(0)  =  0  /  i){x)dx=0.  (3.4) 

J  —  OO 
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a=0,25 

b=3 

J 

k 

1 

Figure  3,1:  The  Gabor  or  Morlet  wavelet 

Therefore,  ip  has  an  average  value  of  0,  acts  like  a  highpass  or  bandpass  filter  and  hence  must  oscillate. 
The  fact  that  ip  e  means  that  ip  has  finite  energy  and  therefore  must  fall  off  fairly  rapidly.  These 

two  features  are  what  motivated  Morlet  and  Grossman  to  refer  to  these  functions  as  “wavelets”  or  short 
waves. 

The  admissibility  condition  (3.1)  was  not  specifically  derived  with  these  features  in  mind.  Rather, 
the  admissibility  condition  is  necessary  for  the  inverse  transform  to  exist  and  for  the  inverse  wavelet 
operator  itself  to  be  a  shifted  and  dilated  version  of  a  mother  wavelet,  referred  to  as  the  dual  of  ip. 

Theorem  2  (Inverse  Wavelet  TVansform)  []0]  Let  ip  be  a  basic  wavelet  which  defines  a  CWT,  Wf. 
Then 

I  roo  yoo  .  .  .dadb 

/(^)  =  7T-  /  /  Wf{a,b)ipa,bix)—Y-  (3.5) 

^—00^—00  d 

□ 


Kaiser  [19]  provides  a  very  understandable  derivation  of  both  the  Inverse  Wavelet  Transform  and  the 
admissibility  condition  that  is  worth  including  here,  since  it  will  be  used  later  to  determine  the  conditions 
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for  dyadic  wavelets  and  frames.  Using  Parseval’s  theorem,  (3.2)  can  be  rewritten  as 

Wfia,  b)  =  (/,  ^a,b)  =  {F,  ^a,b)  (3.6) 

where  F(0  is  the  Fourier  Transform  of  f{x)  and  ^^0,6(0  is  the  Fourier  Transform  of  ^„,6(a;),  given  as 

^a,6(0  =  (3.7) 

where  ^ (^)  is  the  Fourier  Transform  of  ip{x).  Expanding  the  inner  product  of  (3.6)  gives 

/OO  - - - 

-OO 

I  fOO  _ 

=  l«|5  /  (3.8) 

*/ — OO 

The  right  side  of  (3.8)  is  the  Inverse  Fourier  Transform  of  F{()W{^,  so  taking  the  Fourier  Transform 
of  both  sides  with  respect  to  b  gives 

/OO  _ 

Wf{a,b)e~^‘^^^'’db=  |a|2F(^)4r(aO  (3.9) 

'OO 

Kaiser  uses  some  clever  manipulation  to  isolate  F{^)  in  (3.9).  He  multiplies  both  sides  by 
and  integrates  with  respect  to  a  where  the  measure  of  integration  is  da/\a\^. 

J-ooJ-oo  \aY  7-00  |a| 

/OO  rin 

l^(aOPn 

-OO  \Q>\ 

=  FiOZiO  (3.10) 

where 

/OO 

l»K)ln-  (3-11) 

-OO  |a| 

Choosing  da/\a\'^  as  the  measure  associated  with  the  integral  in  (3.10)  guarantees  that  the  admissibility 
condition  is  a  constant  and  therefore  guarantees  that  the  dual  can  be  represented  as  a  shifted  and  dilated 
version  of  a  mother  wavelet  [19].  Now  F{^)  in  (3.10)  can  be  solved  if  and  only  if  0  <  <  Z-i(^)  < 

B  <  00, 

FiO  =  Z-\0  r  r  H^^(a,6)|ap«r(a^)e-i2-«o“.  (3.12) 

^0  J~oo  ar 
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Substituting  u)  —  in  (3. 1 1)  gives  the  admissibility  condition,  which  is  a  constant  and  therefore  bounded. 


Zio  = 


/: 


I^MP 

kl 


duj 


(3.13) 


Substituting  Z(^)  =  in  (3.12)  and  taking  the  Inverse  Fourier  Transform  of  both  sides  with  respect 
to  ^  gives 


I  f<X>  f<X>  /T.-h\ 

/(^)  =  ^/  /  Wf{a,b)\ar2i,(  —  ) 

J~oo  J—oo  \  d  / 


dddl) 


(3.14) 


which  is  the  form  of  the  Inverse  Wavelet  Transform  given  in  Theorem  2. 

Figure  3.2  shows  an  example  of  a  continuous  wavelet  transform.  The  signal  being  analyzed  is  a  tran¬ 
sient  sinusoid  with  exponentially  decaying  amplitude.  The  analyzing  wavelet  is  Morlet’s  wavelet  from 
Figure  3.1.  The  horizontal  axis  is  shift,  6,  and  the  vertical  axis  is  scale,  a  with  a  increasing  downward. 


Figure  3.2:  Continuous  Wavelet  Transform  of  a  transient  signal 


The  transient  signal  has  an  average  value  of  0,  so  for  large  values  of  a,  the  CWT  vanishes  since  a  very 
wide  wavelet  tends  to  average  over  the  transient  signal.  The  CWT  also  vanishes  for  very  small  values  of 
a  since  the  signal  appears  constant  over  a  very  short  interval.  There  is  a  value  of  a,  however,  that  coin¬ 
cides  very  nicely  with  the  oscillation  of  the  transient  signal,  thereby  producing  a  very  large  response  in 
the  CWT,  which  can  be  seen  in  Figure  3.2.  The  time  localization  feature  of  the  CWT,  provided  through 
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the  shift  parameter,  b,  gives  the  approximate  location  of  the  transient  signal. 

The  continuous  wavelet  transform,  characterized  by  a,  b  €  3?,  where  a  ^  0,  is  rather  simple  to  use 
and  requires  very  few  restrictions  on  the  analyzing  wavelet.  However,  it  contains  a  lot  of  redundancy  and 
is  computationally  intensive.  Notice  in  Figure  3.2  the  high  degree  of  correlation  in  the  CWT.  This  cor¬ 
relation  is  due  to  the  redundancy  and  can  be  removed  by  less  redundant  representations.  The  following 
sections  will  show  the  development  of  different  classes  of  wavelets  with  lesser  degrees  of  redundancy 
as  a  and  b  are  further  constrained.  While  these  constraints  add  conditions  on  the  wavelets,  they  lead  to 
the  generation  of  bases,  new  design  techniques,  and  fast  algorithms  for  the  wavelet  transform. 

3.2  Dyadic  Wavelets 

In  order  to  reduce  the  computational  burden  of  the  CWT,  let  the  scale  parameter  take  on  values  of  a  =  2^ 
where  j  ^  Z  and  6  6  3?.  The  wavelet  transform  becomes 

Wfi2^,b)  =  (/,V’2.,6) 

f°°  i  ■  ■ 

=  /  fix)  2  2^(2  ^x-b2  ^)dx  (3.15) 

J -00 

This  class  of  wavelets  is  known  as  “dyadic”  wavelets  and  are  defined  as  follows: 

Definition  4  (Dyadic  Wavelets)  [10]  A  function  ip  €  ij^(3t)  is  called  a  dyadic  wavelet  if  there  exist 
two  positive  constants  A  and  B,  with  0  <  A  <  B  <  oo,  such  that 

OO  2 

E  1^(2''^^) I  (3.16) 

j—-oo 

Condition  (3.16)  is  known  as  the  stability  condition  and  again  is  driven  by  the  necessity  for  an  inverse 
transform.  Since  (3.15)  has  the  same  form  as  the  CWT,  but  with  a  =  2^  then  (3.11)  must  still  hold 
where  a  =  2K  Substituting  Aa  =  2^'^^  —  2^  =  2^  for  da  gives  da/a  =  1.  Summing  over  j  gives 

^(0  =  (3.17) 

jez 

As  before,  Z(^)  must  be  bounded,  which  leads  to  the  stability  condition 

0  <  A  <  k  B  <  oo.  (3.18) 

jez 
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Notice  this  time  that  Z  ((^)  must  not  necessarily  be  a  constant.  However,  since  it  is  bounded,  the  Inverse 
Wavelet  Transform  can  be  found  by  substituting  a  =  li?  and  Aa  =  2^  into  (3.12)  and  summing  over  j 
instead  of  integrating  over  a, 

=  E  (3.19) 

where  =  ^(^)/Z(^)  and  Z(^)  is  given  in  (3.17).  Taking  the  Inverse  Fourier  Transform  of  both 
sides  of  (3 . 1 9)  gives  the  expression  for  the  Inverse  Wavelet  Transform  using  dyadic  wavelets, 

f(^)  =  E  r  2-^Wf(2^,b)2-H  (^)  db  (3.20) 

j^zJ-oo  \  2.J  J 

where  ^(s)  is  the  Inverse  Fourier  Transform  of  #(^)  and  is  the  dual  wavelet  to  V-’(a:). 


3.3  Frames 

The  next  step  in  reducing  the  computational  burden  of  the  CWT  is  to  sample  the  shift  parameter,  letting 

=  k2^bo,  where  bo  is  a  constant  known  as  the  sampling  rate  [10].  Now  the  wavelet  operator  is  given 
as 

^bo;j,k(x)  =  2~2i/;(2-^x  -  kbo)  (3.21) 

and  the  wavelet  transform  given  by 

(aj,  bj,k)  =  if,  'ipbo-,j,k}  (3.22) 

The  condition  on  the  wavelet  also  tightens  beyond  the  stability  condition,  namely, 

^ll/|P<  E  \{f^'^bo-,j,k)f<B\\ff  (3.23) 

j,kez- 

where  ||  •  |p  is  the  L^(3fJ)  norm  and  0  <  .A  <  5  <  oo.  This  condition  is  identical  to  that  for  a  frame  of 
L  (3?),  meaning  that  in  order  for  ^  to  be  a  wavelet  with  the  stated  conditions  on  a  and  b,  it  must  generate 
a  frame  of  L^(3?)  [10]. 
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Definition  5  (Frames  of  L^(3i?))  [10]  A  function  xp  G  is  said  to  generate  a  frame  {V’6o;i,fc} 

L  (3i)  with  sampling  rate  bo  >  0  if  (3.23)  holds  for  some  positive  constants  A  and  B,  which  are  called 
frame  bounds.  If  A  =  B,  then  the  frame  is  called  a  tight  frame. 

Before  giving  the  expression  for  the  inverse  transform,  let  T  be  a  linear  operator  on  defined  by 

'^f  ~  if i'^bo;j,k)'^bo;j,k  (3.24) 

j,keZ 

where  /  G  1/^(3?).  Then  the  dual  of  V’6o;i,fc  is  given  by  and  the  inverse  wavelet 

transform  by 

fi^)=  {f,i’bo-,lk)ip%  (3.25) 

jMZ 

If  A  —  B  =  1,  then  V’6o;i,fc  is  orthonormal  basis  of  Z^(3?).  If  A  =  B  1,  then  the  frame  is  called 
a  tight  frame,  which  acts  like  an  orthonormal  basis,  but  may  not  even  be  linearly  independent  [12].  The 
intent  of  discretizing  both  the  scale  and  shift  parameters  is  to  reduce  the  redundancy  of  the  wavelet  trans¬ 
form.  Frames  provide  an  intermediate  step  between  the  continuous  wavelet  transform,  which  contains 
the  maximum  amount  of  redundancy,  and  wavelet  orthonormal  bases,  which  generate  decompositions 
with  no  redundancy.  The  ratio  of  the  frame  bounds,  B/A,  acts  like  a  redundancy  indicator.  For  instance, 
in  speech  processing,  B/A  is  very  large  indicating  a  lot  of  redundancy  in  the  wavelet  transform  and  in 
fact  approximates  the  continuous  wavelet  transform  [12].  At  the  other  extreme,  applications  like  image 
compression,  requiring  no  redundancy,  constrain  the  wavelets  so  that  they  generate  a  tight  frame,  that 
is,  B/A  =  1. 

A.  Grossman  with  the  help  of  Y.  Meyer  first  realized  the  importance  of  the  frame  concept  with  re¬ 
gards  to  wavelet  analysis.  Meyer  showed  how  the  wavelet  frame  construction  was  the  same  as  that  of 
the  Weyl-Heisenberg  coherent  states.  However,  for  the  W-H  coherent  states,  if  a  basis,  g,  was  required 
(no  redundancy)  as  opposed  to  a  frame,  then  either  xg{x)  or  was  not  square  integrable  (not  in 

L  (3fJ))[12].  Meyer  set  out  to  show  that  Grossman’s  wavelet  frames  were  subject  to  the  same  limitation, 
but  instead  discovered  a  bandlimited,  orthonormal  wavelet  basis  such  that  x'ij){x)  and  were  both 

in  L^(3f?)  [13].  This  discovery  led  many  of  the  wavelet  pioneers  toward  developing  the  mathematics  for 
orthonormal  wavelet  bases. 
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3.4  Orthonormal  Wavelets 


Y.  Meyer  showed  that  for  a  —  2^  and  b  =  k2^  (assume  bo  =  1),  there  exists  a  family  of  functions, 
{^j,k  =  2-Jf'^ip{2-^x  -k)}  that  form  an  orthonormal  basis  of  Chui  and  Mallat  define  an  or- 

thonormal  wavelet  basis  as  follows. 


Definition  6  (Orthonormal  Wavelet  Bases)  [10,  23,  27]  A  function  xp  e  L^i^)  is  called  an  orthonor¬ 
mal  wavelet  if  the  family  is  an  orthonormal  basis  ofL'^{lSi),  that  is, 

{'4’j,ki'^l,m)  —  ^j,l  •  Sk,m  (3.26) 


where 


for  j  =  k 
j^k 


and  every  f  €  1/^(3?)  can  be  written  as 


(3.27) 


00  OO 

/W=  £  E  dl2^{2^x-k)  (3.28) 

j=~oo  k=~oo 

Equation  (3.26)  indicates  that  the  family  of  wavelets,  ^  is  orthogonal  in  two  dimensions.  Integer 
translates  of  the  wavelet  at  a  given  scale  are  orthogonal  to  one  another  (4,^)  as  are  wavelets  at  two 
different  scales  (4,«)-  At  a  given  scale,  j,  {xpj^k,  =  4,m  and  forms  an  orthonormal  basis  of 
Wj,  a  subspace  of  iy2(3{).  Because  of  orthogonality  across  scales,  Wj  1  Wi  for  all  j  ^  1.  Therefore, 
applying  (2.13)  to  the  projection  of  /  onto  Wj  gives 

00 

g^x)  =  {Q!'^f){x)  =  dl2^{2^x  -  k)  (3.29) 

k=—oo 

where  g^{x)  €  Wj  and 

/OO  . _ 

f{x)22xp{2ix-k)dx  (3.30) 

“OO 

is  the  projection  operator  with  respect  to  the  basis  xp.  Substituting  (3.29)  into  (3.28)  gives 


OO 

m  =  E  Kf)(x) 

j=-oo 


(3.31) 
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which  says  that  /  can  be  decomposed  into  an  infinite  sum  of  its  projections  onto  orthogonal  subspaces, 
Wj  [10].  Furthermore,  can  be  expressed  as  the  direct  sum  of  orthogonal  subspaces  Wj 

L2(3i)  =  ...+l^_i+lFo-t-TFi+...  (3.32) 

Because  the  wavelet  has  the  characteristic  of  a  bandpass  filter  (3.4),  the  projection  operator,  is 
effectively  projecting  or  filtering  out  the  detail  or  high  frequency  content  of  /  at  scale  j,  and  (3.3 1 )  shows 
that /consists  of  the  infinite  sum  of  these  detail  functions,  (Qj,/)(x)  [10, 12, 23, 27].  Furthermore,  the 
detail  functions  (Q^/)  (a:),  contained  in  Wj,  are  orthogonal  to  one  another  since  Wj  1  Wi  for  all  j  ^  I 

3.5  Multiresoliition  Analysis 

A  major  breakthrough  in  the  understanding  of  orthonormal  wavelet  bases  came  when  Y.  Meyer  and  S. 
Mallat  imposed  the  concept  of  multiresolution  analysis  on  wavelet  decompositions  [14].  Mallat  had 
been  working  with  the  Laplacian  pyramid  algorithm  developed  by  Burt  and  Adelson  [6]  and  recognized 
that  the  sequence  of  functions  generated  by  an  orthonormal  wavelet  were  in  essence  “detail”  functions 
that  representated  the  information  lost  in  going  from  one  scale  to  a  lower  scale  [13].  He  and  Y.  Meyer 
developed  the  mathematics  for  what  came  to  be  known  as  a  multiresolution  analysis. 

Let  Vj  be  defined  as  a  subspace  of  L^(3?)  where 

Vj  =  ...  +Wj^^+Wj-2\-Wj-i  (3.33) 

Then,  the  direct  sum  decomposition  of  in  (3.32)  can  be  rewritten  [10,  23]  as 

L^(U)  =  Vj+Wj+Wj+i+ . . .  (3.34) 

From  (3.33),  it  is  clear  that 

Vj+I  =  Vj+Wj  (3.35) 

and  Vj  forms  a  nested  sequence  of  subspaces  of  that  is, 

•  •  •  '^J-I  C  Vj  C  Vj+i . . .  (3.36) 
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Let  (t)j^k  be  a  set  of  functions  that  spans  the  subspace  Vj  [10],  where 

(f)j,k  =  22(l){2^x-k)  (3.37) 

Then,  any  function  in  Vj  can  be  represented  by  a  linear  combination  of  (f)j^k- 

OO 

P{x)=  4^^H2^x-k)  (3.38) 

k=~oo 

Assume  <f>j^k  forms  an  orthonormal  basis  of  Vj,  then 

{^j,ki  4>j,m)  =  ^k,m  (3.39) 

and, 

4  =  {fHx),(l}j,k)  (3.40) 

This  assumption  is  not  as  bold  as  it  may  seem.  The  Fourier  transform  of  (3.39),  called  the  Poisson  sum¬ 
mation,  is  a  27r  periodic  function  given  by 

OO 

^  27rfe)p  =  1  (3.41) 

A:=“Oo 

where  $(a;)  is  the  Fourier  transform  of  <f){x).  Let  ^{x)  be  a  function  such  that  spans  Vj,  then  (px  {x) 
can  be  found  such  that  it  satisfies  (3.39)  by  way  of  the  following  orthogonalizing  “trick”  and  taking  the 
inverse  Fourier  transform  [10]: 

- - r  (3.42) 

(Er=-ool^(c^  +  27rfc)P)5 

where  0  <  ^4  <  |  +  27rA:)p  <  .B  <  oo  for  all  a;.  Because  a  non-orthogonal  function 

that  spans  Vj  can  be  orthogonalized  using  (3.42),  it  is  safe  to  assume  that  (l)j^f^  is  orthonormal  in  the  first 
place. 

Now,  if  is  a  function  in  the  subspace  Vj,  then  by  (3.35)  it  can  be  represented  by  the  sum  of  its 
projections  on  Wj-i  and  [10,  12] 

fix)  =  V^i\P){x)  +  Q{^-\P){x)  (3.43) 
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Since  represents  the  detail  of  p  at  scale  j  -  1,  then  must  represent  p  with  those 

details  removed,  that  is,  a  lower  resolution  approximation  of  fj  at  scale  j  -I  [10,  12, 23],  Subsequent 
projections  onto  complement  subspaces,  V  and  W,  produce  the  decomposition  map  shown  in  Figure  3.3. 
Vj  contains  lower  and  lower  resolution  approximations  of  /•?  as  j  -oo.  The  nested  sequence  of 


V.,=L2(9{) 


"j-s 


Wj., 


Figure  3.3:  Multiresolution  Decomposition  of 

subspaces  is  called  a  multiresolution  analysis  and  satisfies  the  following  conditions: 

Definition  7  Multiresolution  Analysis/70,  12,  23,  24,  27] 

The  sequence  of  subspaces,  Vj,form  a  multiresolution  analysis  if  the  following  conditions  are  satisfied: 


1.  . . .  C  F_i  C  Fo  C  Fi . . .  ; 

2.  dosL,  {UjezVj)  =  Lp^); 

=  {0}  >■ 

^+1  =  €  Z;  and 
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5.  f{x)  eVj  ^  f{2x)  e  Vj+i,  j  e  z 


Items  7.1  and  7.2  state  that  a  multiresolution  analysis  is  a  nested  sequence  of  subspaces  whose  union 
spans  1,2(3?).  Item  7.3  states  that  there  is  no  portion  of  1,2(3?)  that  is  common  to  all  subspaces  Vj,  except 
the  all  zero  function,  {0}.  Item  7.4  states  that  a  subspace,  Vj  consists  of  the  sum  of  the  subspace  Vj-i 
and  its  complement  Wj-i.  Finally,  Item  7.5  states  that  if  a  function,  /,  resides  in  the  subspace  at  scale 
j,  then  /  dilated  by  2  resides  in  the  subspace  at  scale  j  +  1,  which  can  be  seen  easily  in  (3.43)  by  letting 
x'  =  2x. 

An  orthonormal  multiresolution  analysis  (MRA)  [10,  12,  23,  24,  27]  further  requires  that  the  two 
complementary  subspaces,  Vj  and  Wj  be  orthogonal  complements,  Vj  ±  Wj,  which  leads  to  the  fol¬ 
lowing  condition  on  their  bases 

=  0  (3.44) 

Let  be  some  function  that  can  be  completely  represented  by  the  orthonormal  basis  of  Vj+i,  that 
is,  /•i+^(a:)  €  Vj+i.  From  (3.43),  can  be  decomposed  into  a  detail  function  and  low  resolution 
approximation  at  scale  j  [10, 12,  23],  that  is. 


=  Vi{f+^)ix)  +  Q^;^{f+^){x)  (3.45) 

=  Pix)+9Hx)  (3.46) 

where  and  reside  in  orthogonal  subspaces,  Vj  and  Wj.  In  a  MRA,  f^{x)  can  be  decomposed  again 
into  orthogonal  subspaces,  Vj-i  and  Wj-i, 


P{x)  =  Vi-\P){x)  +  Q?-\p){x) 

=  P~^{x)  +  p~^{x) 


Substituting  (3.48)  into  (3.46)  gives  a  two  level  decomposition  of 


P^\x)  =  p-\x)+p-\x)^9Px) 
The  A'^-level  decomposition  of  fpx)  is  given  as 

p+\x)  =  P^^-^{x)  +  ^p-^-^\x) 

i=\ 


(3.47) 

(3.48) 


(3.49) 


(3.50) 
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Figure  3.4:  N-level  Multiresolution  Decomposition 

and  is  shown  in  Figure  3.4.  It  is  helpful  to  visualize  what  a  multiresolution  decomposition  looks  like 
by  looking  at  an  example.  Let  the  wavelet  be  Daubechies’  D4  wavelet,  a  known  orthonormal  wavelet, 
shown  in  Figure  3.5  with  its  corresponding  scaling  function.  Daubechies’  technique  for  deriving  these 


V(X) 

_ d 

'T 

Figure  3.5:  Daubechies’  D4  Wavelet  and  Scaling  Function 

wavelets  will  be  discussed  in  more  detail  in  Chapter  4.  Figure  3.6  shows  the  multiresolution  decomposi¬ 
tion  of  a  transient  signal  using  the  D4  wavelet  and  corresponding  scaling  function.  The  functions  on  the 
left  are  the  successive  low  resolution  approximations  of  the  transient  signal  at  each  scale.  The  functions 
on  the  right  are  the  detail  functions  found  by  projecting  the  signal  onto  the  orthonormal  wavelet  basis 
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for  that  scale.  The  basis  is  constructed  by  integer  shifts  of  the  wavelet  dilated  by  2^.  Notice  the  effect 
of  projection.  As  the  signal  is  projected  into  lower  and  lower  scales,  both  the  approximation  and  detail 
functions  begin  to  look  like  the  scaling  function  and  wavelet,  respectively,  which  is  expected  since  each 
is  a  linear  combination  of  the  basis  functions  of  each  subspace,  V  and  W.  This  effect  is  one  reason  why 
one  might  desire  a  wavelet  basis  that  “looks”  like  the  signal  being  analyzed.  The  original  signal  in  Fig¬ 
ure  3.6  is  completely  represented  by  the  detail  functions  and  the  last  low  resolution  approximation  and 
can  be  reconstructed  simply  by  summation. 

As  mentioned  previously,  the  constraints  on  the  wavelet  increase  in  complexity  as  one  moves  from 
the  integral  wavelet  transform  to  orthonormal  bases.  It  is  important  to  understand  those  constraints,  since 
they  will  be  used  extensively  in  the  wavelet  design  algorithm  presented  in  Chapter  5. 


3.5.1  Properties  of  0  and  V' 


As  shown  previously,  in  an  orthonormal  multiresolution,  (j)j^k  is  an  orthonormal  basis  of  the  sub¬ 
space  Vj,  and  is  an  orthonormal  basis  of  the  subspace  Wj,  where  Vj  J.  Wj.  Furthermore,  the  family 
of  functions  {ipj,ki  —oo  <  j  <  oo}  is  an  orthonormal  basis  of  1(^(3?).  These  conditions  lead  to  the  fol¬ 
lowing  relationships  between  (f)j^k  and  [10] 


II 

(3.51) 

o 

II 

(3.52) 

(3.53) 

as  given  previously  in  (3.39)  (3.44)  and  (3.26).  Furthermore,  since  ip{x)  is  a  wavelet,  then  by  (3.4) 


/: 


'>p{x)dx  =  0 


'J'(O)  =  0 


(3.54) 


By  (3.35)  it  is  clear  that  since  <f)j^k  and  'ipj^k  are  bases  of  Vj  G  l^+i  and  Wj  G  T^+i,  respectively, 
they  both  reside  in  and  can  therefore  be  represented  by  a  linear  combination  of  the  basis  of 
[10, 12, 23, 31].  For  j  =  0,  (/>o,o  =  G  Vj  can  be  represented  using  (3.38)  as: 

OO 

((){x)=  Y,  Ck2^(j){2x-k)  (3.55) 

k=-oo 
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Since  V’(a;)  also  resides  in  Vi,  (3.38)  gives  the  generating  equation  for  the  wavelet 

OO 

i’{x)=  h^hi^x-k).  (3.56) 

k=—<X) 

The  2^/2  term  in  both  of  the  above  equations  can  be  included  in  the  coefficients  giving 

OO 

X]  Pk(l>{2x-k)  (3.57) 

k=~oo 

and 

OO 

^{x)=  QkH'^x-k)  (3.58) 

fc=  — OO 

Equation  (3.57)  is  a  recursive  equation  called  the  two  scale  relation  for  (f>,  and  pk  and  qk  are  called  the 
generating  sequences  for  (j)  and  'ip,  respectively  [10].  Taking  the  Fourier  Transform  of  (3.57)  and  (3.58) 
gives  the  relationship  between  $(c4;),  ^'(a;),  P(c<;),  and 


(3.59) 


(3.60) 


Assume  pix),  called  the  scaling  function  [10,  12,  23,  24,  27],  generates  the  multiresolution  analysis, 
{Vj},  and  is  normalized  [31]  such  that 


/OO 

(l>{x)d 

-OO 


lx  =  1 


<&(0)  =  1 


(3.61) 


Since  ^(x  +  n)  and  'ip{x  +  n)  form  orthonormal  bases  of  Vq  and  Wq,  respectively,  and  are  orthogonal 
to  one  another,  then  the  following  must  be  true 


fOO 

/  \(j){x)\^  dx  =  1 

J  —  OO 

(3.62) 

fOO 

/  \'ip{x)f  dx  =  1 

J—oo 

(3.63) 

^{x)(j){x  +  n)dx  =  5{n) 

0 

(3.64) 

‘ip{x)tp(x  +  n)dx  =  S{n) 

D 

(3.65) 

OO 

(f){x)'ip{x  +  n)dx  =  0 

-OO 

(3.66) 
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Taking  the  Fourier  Transform  of  (3.64)  gives  the  Poisson  Summation  on  $(cj)  [23,  24] 


T 


^{11  (t>{x)(f>{x  +  n)dx^ 

<  00  \ 

[(j){x) -k  (f){x)]  ■  ^  5(a;  +  n)| 

k  n=— oo  / 

oo 

|$(a;)p  *  ^  6{u!  +  2Trn) 


=  nS{n)) 
=  1 

=  1 


7l=— 00 
00 


y~!  |$(u;  +  27rn)p  =  1 


(3.67) 


n——oo 

where  is  the  correlation  operator  and  is  the  convolution  operator.  Likewise,  taking  the  Fourier 
Transform  of  (3.65)  gives  the  Poisson  Summation  on  5f(u;). 


OO 

|4'(a;  +  27rn)|^  =  1 

n=“00 


(3.68) 


Substituting  (3.57)  and  (3.58)  into  (3.4)  and  (3.61)-(3.66)  produces  several  conditions  on  pk  and  [10, 


12] 


Two  additional  conditions  on 
are  provided  in  Section  3.6. 


°°  nr 

Pk  =  2  ^  P(0)  =  2 

(3.69) 

k=—oo 

OO  _ 

^k  =  0  Q{0)  =  0 

(3.70) 

k=-oo 

00 

E  pI  =  2 

(3.71) 

k—~oo 

oo 

*ik  =  ^ 

(3.72) 

k=—oo 

00 

E  PkPk-2n  =  26 {n) 

(3.73) 

k~—oo 

00 

E  *lkQk-2n  =  26{n) 

(3.74) 

/c  =  — 00 

oo 

E  PkQk-2n  =  0. 

(3.75) 

k=—oo 

P{uj)  and  Q{uj)  are  stated  here  without  proof  or  explanation.  The  details 

l^(w)P  +  |Q(w)p  =  4  (3.76) 
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P{u})P{ijj  +  tt)  +  Q{(jj)Q{ij0  4-  tt)  =  0 


(3.77) 


Conditions  (3.76)  and  (3.77)  allow  for  perfect  reconstruction  of  the  decomposition  of  /. 

It  is  important  to  note  that  the  scaling  function,  (j){x),  generates  the  multiresolution  analysis,  {V}} 
and  the  wavelet,  'ij){x),  generates  the  multiresolution  decomposition,  {g^  (a:)}.  Before  the  MRA  was  for¬ 
mulated  by  Meyer  and  Mallat,  there  was  no  standard  technique  for  finding  orthonormal  wavelet  bases. 
Through  the  scaling  function,  however,  several  design  techniques  were  developed  and  they  will  be  dis¬ 
cussed  in  Chapter  4.  Furthermore,  because  the  scaling  function  is  the  MRA  generator,  the  conditions  for 
an  orthonormal  MRA  rest  on  the  scaling  function  as  will  be  shown  in  Chapter  5. 

3.6  Discrete  Wavelet  TVansform 

The  purpose  for  moving  from  the  continuous  wavelet  transform  (CWT)  to  frames  and  ultimately  to  or¬ 
thonormal  bases  was  to  remove  the  redundancy  in  the  wavelet  transform.  Multiresolution  analyses  sum¬ 
marized  in  Section  3.5  provide  the  framework  for  finding  orthonormal  bases  of  L^(3J),  but  are  still  based 
on  continuous  functions.  That  is,  the  orthogonal  multiresolution  decomposition  given  in  (3.50)  starts 
with  a  continuous  function  and  decomposes  it  into  a  series  of  detail  functions  and  a  final  low  resolution 
residual.  However,  each  detail  and  low  resolution  function  can  be  uniquely  defined  by  its  correspond¬ 
ing  projection  coefficients,  4  and  4  at  some  scale  j,  as  shown  in  (3.29),  (3.38)  and  Figure  3.4.  In  some 
digital  signal  and  image  processing  applications,  like  compression,  one  would  like  to  deal  with  the  co¬ 
efficients  only  and  never  have  to  produce  the  continuous  detail  or  low  resolution  functions.  Mallat  [24] 
used  his  MRA  construction  to  develop  a  fast  algorithm  for  finding  the  multiresolution  decomposition  of 
a  signal.  In  Figure  3.4,  let  =  /,  then 

oo 

k=~oo 

oo 

~  (3.78) 

k=—oo 
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and  likewise, 


k=~oo 


(3.79) 


Through  a  change  in  variables  and  some  algebra  it  can  be  shown  [10, 12,  23]  that 

-i  1 

=  22/  <p{-x)(f){x  -  {k-  2m))dx 

J~oo  ^ 


(3.80) 


Notice  that  the  inner  product  of  two  scaling  functions  from  adjacent  scales  is  not  a  function  of  their  scales 
at  all!  The  scale  parameter,  j,  does  not  appear  in  the  right  side  of  (3.80).  Similarly, 


Let 


and 


=  2“ 2  /  x)<f){x  —  {k  —  2m))dx 

J~oo  ^ 

/OO 

(l){-x)(f>{x  —  n)dx 

-OO  ^ 


then  (3.78)  and  (3.79)  become 


and 


-1  f°°  1 

9n  =  ^  /  i’{-x)(j){x  -  n)dx 

J  —  OO  ^ 


^  ^  ^k^k—2m 
k——oo 


=  Yi  49k-2m 

k——oo 


(3.81) 


(3.82) 


(3.83) 


(3.84) 


(3.85) 


where  the  remaining  2  term  is  included  in  c^.  Equations  (3.84)  and  (3.85)  show  that  the  projection 
coefficients  for  the  low  resolution  and  detail  functions  can  be  obtained  by  filtering  and  downsampling 
the  low  resolution  projection  coefficients  from  the  previous  scale.  Furthermore,  the  digital  filters  used 
are  the  same  regardless  of  which  scale  is  being  decomposed.  Equations  (3.84)  and  (3.85)  are  known  as 
the  Discrete  Wavelet  Transform.  [10,  12,  23,  30]  The  multiresolution  decomposition  shown  in  Figure 
3.4  can  be  represented  by  a  sequence  of  Discrete  Wavelet  Transforms  (Figure  3.7)  [23].  Substituting 
(3.57)  and  (3.58)  into  (3.82)  and  (3.83),  respectively,  gives 


K  =  2^71 


(3.86) 
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Figure  3.7:  Discrete  Wavelet  Transform 
1 

9n  =  -j^qn-  (3.87) 

So,  except  for  a  constant  scale  factor,  the  digital  filters,  hn  and  used  to  implement  the  Discrete  Wavelet 
Transform  are  precisely  the  generating  sequences  for  <f){x)  and  ip{x),  respectively  [12].  Given  this  rela¬ 
tionship  and  the  fact  that  for  digital  signals,  one  would  like  to  deal  with  only  the  projection  coefficients 
and  the  digital  filters,  as  in  (3.84)  and  (3.85),  the  conditions  on  p  and  q  in  (3.59),  (3.60),  (3.69),  (3.70), 
and  (3.71)-(3.77)  can  be  translated  directly  to  conditions  on  h  and  g. 


OO 


=  2  ^  hk(l){2x  -  k) 

k=—oo 

(3.88) 

OO 

ip{x)  =  2  Y,  9k<t>{2x-k) 

k=—oo 

(3.89) 

♦M  =  (1) «  {f) 

(3.90) 

»(")  =  (^)  ♦  (1) 

(3.91) 

OO 

E  hk  =  l 

k=—oo 

(3.92) 

o 

II 

(3.93) 

k=—oo 
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(3.94) 


00  ^ 

E  >‘1  =  1 

k=~oo 

oo  ^ 

E  9*  =  5 


hkhk~2n  —  2^(^) 

k~—oo 


9k9k—2n  — 

A;=-oo  ^ 


hk9k~2n  —  0 

k=—oo 

\H{u)'^  +  \G{uj)\‘^  =  1 


H{u;)H{oj  +  tt)  +  G{u)G{uj  +  tt)  =  0 


(3.95) 

(3.96) 

(3.97) 

(3.98) 

(3.99) 
(3.100) 


From  (3.92)  and  (3.93)  it  can  be  shown  that  hk  and  gk  are  low  pass  and  high  pass  filters,  respectively 
[10,  12,  23,  30].  From  (3.84)  and  (3.85),  is  found  by  passing  through  a  low  pass  filter  and 
by  passing  through  a  high  pass  filter,  which  is  consistent  since  represents  a  low  resolution  ap¬ 
proximation  of  the  original  sequence  and  represents  the  details  (or  high  frequency  content)  of  the 
original  signal. 

Since  (3.90)  is  recursive,  (3.90)  and  (3.91)  can  be  rewritten  [10,  12,  23,  24]  as 

OO  y  V 

(3.101) 

*M  =  G(^)n«(^)  (3.ro2) 

Continuing  with  Daubechies’  D4  wavelet  as  an  example.  Figure  3.8  shows  the  same  multiresolution 
decomposition  as  Figure  3.6  but  in  terms  of  the  projection  coefficients,  cj.  and  d^,  only.  Notice  that  as  the 
decomposition  proceeds  to  more  levels,  the  number  of  coefficients  decreases  due  to  the  downsampling 
in  (3.84)  and  (3.85). 

In  applications  like  image  compression  or  image  enhancement,  it  is  necessary  to  reconstruct  the  orig¬ 
inal  coefficients  from  the  decomposition  coefficients.  Equations  (3.99)  and  (3. 1 00)  are  the  conditions  on 
H  and  G  that  guarantee  perfect  reconstruction,  and  they  can  be  derived  by  forming  the  reconstruction 
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expression  for4  in  Figure  3.4  [12],  f  is  the  projection  of  /  onto  Vj,  where  /  =  /•?+!  in  Figure  3.4, 
and  are  the  projection  coefficients  of  f  €  Vj, 

4  =  (3.103) 

Substituting  (3.43)  into  (3.103)  where  =  V^f  gives 

OO  oo 

~  ^  4n^{i>j-l,m,4'j,k)  +  ^  d^m^{i’j-l,m,<l>j,k) 

m=— OO  m=— 00 

00  oo 

~  X/  ^  hk-2m+  ^  di^^gk-2m  (3.104) 

m=— oo  m=— oo 

So,  reconstruction  is  accomplished  by  upsampling  the  two  sets  of  coefficients,  and  and  in¬ 
terpolating  with  the  same  filters,  hjn  and  gm,  respectively.  A  single  stage  decomposition/reconstruction 
cycle  is  shown  in  Figure  3.9  [23].  Figure  3.10  shows  the  signal  spectrum  at  each  stage  of  the  decom- 


Figure  3.9:  Decomposition/Reconstruction  cycle  with  QMF  filters 
position/reconstruction  process.  Let  the  Fourier  Transforms  of  the  sequences  in  Figure  3.9  be  given  by 

h4)  =  c^i^) 
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Nyquist 


^ ;  f)  Spectrum  of  and  ;  g)  Spectrum  of 
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In  order  to  have  perfect  reconstruction,  C^u)  must  equal  An  expression  for  C{u;)  is  ob¬ 

tained  by  working  backwards  through  Figure  3.9. 

&{ui)  =  +  D^~\u})G{u;)  (3.105) 

Upsampling  and  to  obtain  ^  and  is  done  by  placing  zeroes  between  each  sample, 
causing  a  simple  dilation  of  the  respective  frequency  spectra  (Figures  3.10d  and  3.10e),  that  is, 

=  C^-^{2uj) 

=  D^-^{2u)  (3.106) 

Downsampling  ^  and  ^  to  obtain  and  is  not  as  straightforward,  because  it  causes  alias¬ 
ing  (Figures  3.10c  and  3.10d)  [35].  The  filtered  spectra,  C^-^u)  and  ^^-^(a;)  are  27r-periodic  and 
oversampled  (Figure  3.10c).  However,  since  H{u)  and  G{lo)  are  not  ideal  filters  (Figure  3.10b),  the 
oversampling  of  &  ^(w)  and  is  less  than  2,  and  therefore,  downsampling  by  2  causes  alias¬ 

ing [35].  The  expression  for  C-?-i(w)  andD^-^u)  intermsof  C-?-i(a;)  andD^-^uj)  must  include  the 
aliasing  terms 

Finally,  C^~^{uj)  and  .D-’“^(a;)  are  filtered  versions  of  the  input  spectrum,  C-^  (a;) 

=  C^{u)H{u) 

D^-Hu)  =  CHu;)G{^ 

The  complex  conjugates  of  the  filter  spectra  in  (3. 108)  are  due  to  the  digital  filters  being  index-reversed, 
h-k  and  g^k-  Substituting  from  (3.108)  back  to  (3.105)  gives 

C[io)  =  ^iy(u;)Ff(a;)  -1-  G(ca)G(w)j  C^{uj)  -\- 

^H{u)H{lu  -1-  tt)  +  G{u)Giuj  +  tt)]  C^{uj  +  tt)  (3.109) 


(3.107) 


(3.108) 
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The  two  conditions  for  perfect  reconstruction  emerge  from  (3.109)  [12,  35] 

|ff(a;)|2  +  |(?(u;)|2  =  l  (3.110) 

H{u)H{u} +  Tr)  +  +  n)  =  0  (3.111) 

The  formulation  of  these  filter  requirements  for  perfect  reconstruction  is  well  known  in  the  field  of  sub¬ 
band  coding  [35]. 

If  the  relationship  between  H{co)  and  G{u)  were  constrained  such  that 

G(u>)  =  e-^H{u  -I-  tt)  (3.112) 

then  the  second  condition  (3.100)  is  always  satisfied  [12,  35]  and  (3.99)  becomes 

\H{u;)f +  \H{u  +  Tr)f  (3.113) 

These  filters  are  called  “quadrature  mirror  filters  (QMF)”  and  the  structure  of  Figure  3.9  is  known  as  a 
2-band  QMF  [35].  Taking  the  inverse  Fourier  Transform  of  (3.112)  gives  the  relationship  between  h  and 
9 

(3.114) 

Substituting  (3.86)  and  (3.87)  into  (3.1 14)  gives  the  relationship  between  the  generating  sequences  that 
guarantees  perfect  reconstruction  in  a  multiresolution  analysis 

(3.115) 

Assume pfc,  a  generating  sequence  for  (j){x),  exists  such  that  (3.69)  and  (3.73)  are  satisfied.  Then,  clearly, 
conditions  (3.92)  and  (3.96)  on  will  also  be  satisfied.  These  conditions  imply  that  is  3n  orthonor¬ 
mal  basis  of  Vj  and  that  /  (l){x)dx  =  1.  Assuming  the  relationship  between  hk  and  in  (3.1 14),  it  can 
be  shown  by  substitution  of  (3.114)  that  every  condition  (3.90)-(3.100)  is  satisfied  and  (j){x)  generates 
an  orthonormal  multiresolution  analysis!  The  two-scale  relation  for  ip{x)  becomes 

00 

'4>{x)=  {-lf^^pi^k4>{‘^x-k)  (3.116) 

k=—oo 
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3.7  2-Dimensional  DWT 


A  2-diniensional  multiresolution  analysis  (MRA)  is  defined  as  a  sequence  of  subspaces,  {Vf},  which 
satisfies  the  conditions  in  Definition  7,  but  defined  on  the  Hilbert  space  [12],  Assume 

-  k^2^y  -  1)  is  an  orthonormal  basis  of  Vf  C  where  (j>{x,y)  is 

the  2-dimensional  scaling  function  that  generates  the  orthonormal  MRA,  {Vf}.  Meyer  showed  that  if 
y)  is  separable,  then  Vf  is  the  tensor  product  of  two  identical  subspaces  of  [12] 

vf  =  V,  ®  Vj 

where  Vj  C  The  2-dimensional,  separable  scaling  function  becomes 

<t>ix,y)  =  Hx)  •  (t>iy)  (3.117) 

where  (pix)  and  ip{y)  generate  identical  orthonormal  multiresolution  analyses,  {Vj},  of  and  4)j,k{x) 
and  0jy(y)  form  orthonormal  bases  of  Fj  c  L'^{U).  Just  as  in  the  1 -dimensional  case,  a  signal, /•?+!  (a;,  2/), 
at  scale  j  +  1  can  be  represented  by  the  sum  of  its  projections  into  orthogonal  subspaces  at  scale  j. 
Given  separability,  P+^{x,y)  is  projected  in  x  onto  Vj  and  Wj  and  then  in  y  onto  Vj  and  Wj.  Since 

~  ^+1  ^+1’  and  Vjj^i  =  Vj+Wj,  then  the  2-dimensional,  separable  projection  produces  4 

subspaces  [12]. 

=  {Vj+Wj)®{Vj+Wj) 

=  Vj  ®  Vj+Vj  ®  Wj+Wj  ®  Vj+Wj  ®  Wj 

where  the  first  subspace  is  V^  and  contains  the  low  resolution  approximation,  P{x,  y),  and  the  remain¬ 
ing  three  subspaces  are  wavelet  subspaces  that  contain  some  version  of  the  details  projected  from 
(x,  y).  The  2-dimensional  bases  for  each  of  these  subspaces  is  given  as  follows  [12]: 

^hk{x)(t)j,i{y) 

<f>j,kix)^j,iiy) 
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^hk{x)i>j,i{y). 


The  multiresolution  decomposition  of  a  2-dimensional  signal  is  shown  in  Figure  3.11.  From  Figure  3.11, 


Figure  3.11:  2-D  Multiresolution  Decomposition 


00  OO 

l=~oo  k=—oo 

OO  OO 


oo  oo 

/=— oo  k=—oo 
00  oo 

=  H  H  4,lhk-2mhl-2n  (3.118) 

/=  — oo  k=—oo 

Likewise, 

oo  00 

^^771,71  ”  X^  ^k,l^k-2m9l~27i  (3.119) 

l=—ook=—oo 
00  oo 

”  X^  X^  ^k,l9k-2mhl~2n  (3.120) 

;=— 00  k~—oo 
00  oo 

“  X^  X^  ^k,l9k~2m9l-27i  (3.121) 

/z=— oo  k~—oo 

Equations  (3.118)-(3.121)  constitute  the  2-dimensional  Discrete  Wavelet  Transform  [12,  23].  The  low 
resolution  approximation  of  the  original  sequence  is  found  by  low  pass  filtering  in  both  the  row  and 
column  dimensions,  is  found  by  low  pass  filtering  the  rows  and  high  pass  filtering  the  columns, 
high  pass  filtering  the  rows  and  low  pass  filtering  the  columns,  and  by  high  pass  filtering 
both  the  rows  and  columns.  Figure  3.12  shows  an  original  256x256  image  and  its  4-level  multiresolution 
decomposition  computed  with  the  2-D  DWT  based  on  Daubechies’  DA  wavelet. 


Figure  3.12:  Multiresolution  Decomposition  of  Lena  Using  Daubechies’  D4  Wavelet 
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3.8  Limitations  to  the  MRA  and  DWT 


3.8.1  MRAs 

While  the  MRA  construction  developed  by  Mallat  and  Meyer  [24]  provides  both  time  and  frequency  lo¬ 
calization  of  signals  containing  both  high  and  low  frequencies,  something  the  Short  Time  Fourier  Trans¬ 
form  (STFT)  was  unable  to  do,  it  does  not  do  a  good  job  of  representing  all  functions  in  1,^(3?)  [25], 
Signals  with  narrowband,  high  frequency  components  are  not  well  represented  because  of  the  constant 
Q  feature  of  the  passbands.  When  the  scale  parameter,  j  in  -  k),  increases  by  1,  the  scale 

doubles  as  does  the  bandwidth.  The  constant  Q  restriction  does  not  provide  independent  control  over 
center  frequency  of  a  passband  and  its  bandwidth.  Mallat  and  Zheng  make  provisions  for  independent 
control  by  including  a  phase  modulation  term  to  their  time-frequency  atom  [25],  Wickerhauser  gener¬ 
alized  Mallat’s  MRA  when  they  developed  the  wavelet  packet  paradigm  [36].  In  an  MRA,  a  signal  is 
projected  into  two  orthogonal  subspaces,  Vj  and  Wj.  Subsequent  decomposition  is  done  on  («)  e  Vj. 
Wickerhauser  removed  this  constraint  and  allowed  the  detail  signal,  g^x)  €  Wj  to  be  decomposed  if  it 
had  the  dominant  energy.  The  resultant  decomposition  tree,  an  example  of  which  is  given  in  Figure  3.13, 
can  take  on  0(2^  —  1)  different  configurations  where  N  is  the  number  of  levels  of  the  decomposition. 
Mallat’s  MRA  is  only  one  of  those  configurations. 

For  many  naturally  occuring  signals,  however,  the  bandwidth  of  a  signal  component  is  a  function  of 
its  center  frequency.  For  example,  a  radar  transmitter  can  transmit  a  signal  with  a  bandwidth  on  the  order 
of  10%  of  its  center  frequency,  in  which  case,  the  Q  of  the  transmitter  would  be  0. 10.  Before  applying 
wavelets  to  an  application,  it  is  important  to  determine  which  class  of  wavelet  is  appropriate. 

3.8.2  DWT 

The  discrete  wavelet  transform  implements  Mallat’s  multiresolution  with  digital  filters.  The  decompo¬ 
sition  equations  for  the  low  resolution  and  detail  projection  coefficients,  and  d{  are  given  as 

OO 

4=  E  4il^hm-2k  (3.122) 

m=— OO 
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V.=L2(3!) 


Figure  3.13:  Example  Wavelet  Packet  Decomposition 


(3.123) 


cx:) 

4=  II  4i'^gm-2k 

m=— 00 

where  h  and  g  are  low  and  high  pass  filters,  respectively.  Notice,  however,  that  neither  (3.122)  nor 
(3.123)  are  shift  invariant.  Let  the  input  sequence  be  shifted  by  n  e  then  (3.122)  becomes 

OO 

4  =  E  4;t\hm-2k 

m=-oo 

OO 

=  E  44^ ^e-2(k-n/2) 
i—-00 

+  4-n  (3.124) 

and  likewise  for  (3.123).  The  fact  that  the  DWT  is  shift  variant  makes  it  very  difficult  to  use  in  applica¬ 
tions  like  object  detection  and  pattern  recognition  because  the  output  of  the  decomposition  is  dependent 
on  the  input  sequence  [23]. 

3.8.3  2-D  DWT 

The  same  limitations  on  the  DWT  apply  to  the  2-D  DWT.  The  multiresolution  decomposition  coeffi¬ 
cients  of  an  image  are  dependent  on  the  input  image,  making  repeatability  for  a  given  class  of  images 
impossible.  Furthermore,  the  implementation  of  the  2-D  DWT  developed  in  Section  3.7  assumes  sepa¬ 
rable  2-D  wavelets  and  scaling  functions.  This  assumption  greatly  simplifies  the  2-D  DWT,  but  results 
in  an  algorithm  that  is  very  sensitive  to  image  rotation.  This  limitation  makes  the  2-D  DWT  even  less 
attractive  for  pattern  recognition  or  object  detection. 
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Chapter  4 


Wavelet  Design  Techniques 


The  development  of  wavelet  multiresolution  analyses  provides  the  necessary  framework  for  designing 
whole  new  classes  of  wavelets.  Ingrid  Daubechies  provided  remarkable  insight  into  wavelet  design 
when  she  published  her  technique  for  finding  orthonormal  wavelet  bases  with  compact  support  [12], 
Others  used  the  mathematical  construct  provided  by  an  MRA  to  define  scaling  functions  and  their  corre¬ 
sponding  wavelets  that  were  optimum  in  some  sense.  Splines  were  used  to  achieve  maximum  smooth¬ 
ness,  projections  were  used  to  find  wavelets  that  were  somewhat  matched  to  the  signal  of  interest.  In 
each  of  these  cases,  however,  the  design  emphasis  is  on  the  scaling  function,  (f>{x)  or  its  generating  se¬ 
quence,  hk,  which  is  to  be  expected  since  the  MRA  is  generated  by  the  scaling  function.  In  this  chapter, 
6  different  design  techniques  will  be  reviewed  in  order  to  provide  motivation  for  the  new  wavelet  design 
technique  presented  in  Chapter  5. 

4.1  Compactly  Supported  Wavelets 

In  1987,  Ingrid  Daubechies,  inspired  by  Mallat’s  Discrete  Wavelet  Transform  algorithm,  succeeded  in 
constructing  orthonormal  wavelet  bases  with  compact  support  in  the  space  domain  [14].  Her  technique 
takes  advantage  of  the  orthonormality  constraints  put  on  the  digital  quadrature  mirror  filter  (QMF),  h^, 
which  is  also  the  generating  sequence  for  the  scaling  function,  <j){x),  of  an  MRA.  It  is  easy  to  see  from 
(5.17)  and  (5.1 8)  that  a  finite  generating  sequence  giyes  rise  to  a  compactly  supported  scaling  function 
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and  wavelet.  Assuming  hk  is  finite  in  length,  M  =  2N  —  1,  its  Fourier  Transform  becomes  a  polynomial 
in 

2A^-l 

H{u:)  =  E  (4.1) 

k=0 

Daubechies  added  constraints  on  H  (a;)  that  would  guarantee  some  degree  of  smoothness  in  the  result, 
namely  that  H (cj)  should  have  N  zeros  at  u;  =  tt,  giving  H  (u;)  the  following  form, 

W{u).  (4.2) 

Since  H{u)  is  a  polynomial  of  degree  2N-1,  Daubechies  is  essentially  dedicating  N  degrees  of  freedom 
in  H{uj)  toward  guaranteeing  smoothness.  The  remaining  task  is  to  find  polynomials  of  the  form  (4.2) 
that  satisfy  the  orthonormality  condition  (3.113).  The  set  of  all  such  polynomials  can  be  found  using 
some  rather  sophisticated  theorems  of  algebra.  Kaiser  [19],  however,  follows  a  shortcut  showing  that 
the  solution  for  iV  =  1,  which  is  the  Haar  system,  actually  generates  solutions  for  all  N  >  2.  For 
N  =l,hk  =  {1/2, 1/2}  giving 


For  N>2,  let 


1  4- 

H{u)  =  =  e*2  cos(a;/2) 

1  — 

H{u  +  tt)  =  — - —  =  sm(w/2). 


C(a;)  =  cos(u;/2)  and  5(0;)  =  sin(a;/2). 


(4.3) 

(4.4) 


(4.5) 


Raising  +  5^  =  1  to  the  power  2N  —  1  gives 


(4.6) 


where  (^)  =  M\/k\{M  —  k)\  are  the  binomial  coefficients.  Kaiser  shows  that  if  the  polynomial, 
?^Ar(C),  consists  of  the  first  N  terms  of  (4.6)  and  H{u!)  can  be  found  such  that  \H{(x>)\'^  =  HNiC), 
then  \H{uj  +  is  automatically  the  last  N  terms  of  (4.6),  thereby  satisfying  the  orthonormality  con¬ 
dition  in  (3.113)  [19].  From  (4.6),  V.  can  be  factored  as  follows. 


n  =  =  C^^Wn{S)  = 


l-t-e*‘^ 

2 


2N 

W;v(5) 


(4.7) 
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where 


^  /2iV-  f 
k  , 


C 


•2N-2-2k  c2k 


(4.8) 


N-1  ( 

Wn{S)  =  E 

fc=0  \ 

The  expression  in  (4.7)  has  exactly  the  form  of  the  squared  magnitude  of  (4.2).  So  all  that  is  left  is 
finding  a  polynomial,  W{u)),  such  that  |VF(a;)p  :=  Wi\f{S{u)),  which  can  be  done  through  algebraic 
manipulation. 

An  example  for  iV  =  2  will  show  the  complete  process, 


=  e(* 

')  (1  -  52)1-'=52'= 

(4.10) 

7 

=  (1-5^)  + 352 

(4.11) 

=  1  +  252 

(4.12) 

=  1  +  2| 

(4.13) 

K  ^-2  ) 

^  g-io.  j 

(4.14) 

Since  W2('5)  =  \W (a;)p  and  is  a  polynomial  in  of  order  2,  then 
2-^(e*‘^  +  e-*‘")  =  |a  +  6e*‘^|2 

=  {a^  +  h'^)  +  ah{e^'^  + 
which  gives  a  =  5(1  +  ■v/S)  and  6  =  ^(1  —  \/3)-  Expanding  (4.2)  gives 

/l  +  eia>\ 


H[u;)  = 


(a  +  be'^) 


i  (fl  +  (2a  +  h)e^  +  {a  +  2by'^  +  6e*3‘^) 


so  the  4  coefficients  for  the  QMF  filter  of  Daubechies’  D4  system  are  given  by 


ho 

a 

1  +  \/3 

hi 

2(1  -|-  b 

_  1 

3  +  y/3 

h2 

a  +  26 

~  8 

3-x/3 

.  '''3  . 

^  6 

1  -  Vs 

(4.15) 


(4.16) 


(4.17) 
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Once  the  filter  coefficients,  hk,  are  found,  ^(x)  can  be  found  recursively  using  (5.17)  and  '4j{x)  can  be 
found  directly  using  (5.18)  where 

This  technique  is  very  clever  and  was  considered  landmark  in  the  development  of  orthonormal  wavelet 

bases.  However,  as  can  be  seen  from  this  entire  process,  the  design  is  focused  on  the  low  pass  filter,  hk- 

> 

Nowhere  in  this  technique  are  the  requirements  on  the  wavelet  addressed.  After  finding  the  generating 
sequence,  h^,  the  scaling  function  and  wavelet  are  found  using  (5.17)  and  (5.18).  The  wavelet  is  driven 
entirely  by  hk,  and  subsequently,  <f){x). 


4.2  Wavelets  for  Signal  Representation 


Tewfik,  Sinha,  and  Jorgensen  [34]  developed  a  technique  for  finding  an  orthonormal  wavelet  with  com¬ 
pact  support  that  provides  the  “best”  signal  representation  of  a  specified  signal  over  a  finite  number  of 
scales.  In  an  orthonormal  multiresolution  analysis  (MRA),  a  signal,  f{x)  £  can  be  represented 
by  an  infinite  sum  of  detail  functions,  (x),  given  in  (3.28)  and  (3.29),  where  the  projection  operator  is 
the  inner  product  with  the  orthonormal  wavelet  basis,  2^/‘^'ip{2^x  -  k).  Truncating  the  infinite  sum  to  N 
scales,  say  scale  jtoj  +  l  —  N  gives  the  following  approximation  of  f{x), 

fix)  =  f+\x)  =  ^pi+i-*(x)  -h  f+^-^ix).  (4.18) 

i=l 

However,  since  —  Wj+Wj-i+ . . .  +Wj-N^i+Vj  -N+1,  (4.18)  is  equivalent  to  /(x)  being  rep¬ 
resented  by  the  orthonormal  basis  of  I'j+i,  namely  —  k),  where  (f){x)  is  the  scaling 

function.  So,  (4.18)  can  be  rewritten  as 

OO 

/(^)=  afe2^(^(2^+ix  -  A:)  (4.19) 

k=—oo 

which  in  the  frequency  domain  is  given  as 

F{lj)  =2~iA{2~^Lj)^{2~^u)  (4.20) 


where  A{uj)  is  the  Fourier  transform  of  ak  and  therefore  27r  periodic.  The  error  norm  between  the 
signal  and  its  approximation  is  given  in  the  frequency  domain  as 

||f(u;)  -  2-2  A(2-^u;)$(2-^w)||  =  ||22F(2^u;)  -  A(a;)$(cu)|| .  (4.21) 
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Tewfik,  et.  al.,  assuming  compactly  supported  wavelets,  parameterize  the  generating  sequence,  h^,  us¬ 
ing  Daubechies  polynomial  factoring  “trick”  and  develop  an  upper  bound  to  the  error  in  (4.21)  where 
the  bound  is  in  terms  of  hk  only.  Once  the  upper  bound  is  given,  they  use  constrained  optimization  tech¬ 
niques  to  solve  for  the  finite  length  hk  that  minimizes  the  upper  error  bound  on  the  signal  representation. 
Once  again,  as  in  Daubechies’  technique,  the  design  emphasis  is  on  hk,  the  generating  sequence  of  the 
scaling  function.  No  design  criteria  are  put  on  the  wavelet,  whatsoever.  This  comes  about  because  as 
shown  in  (4.19),  a  finite  wavelet  decomposition  can  be  recast  as  a  projection  of  the  signal  into  the  next 

higher  approximation  subspace,  Vj^i ,  and  the  only  contributor  to  the  approximation  accuracy  is  the  scal¬ 
ing  function. 

Gopinath,  Odegard  and  Burras  [18]  expanded  the  work  of  Tewfik,  et.  al.,  by  assuming  bandlimited 
wavelets.  This  additional  condition  simplified  the  upper  bound  of  the  error  norm  (4.21)  giving  way  to  a 
much  simpler  numerical  solution.  However,  Gopinath,  et.  al.,  likewise  acknowledge  the  importance  of 
the  scaling  function  in  their  approach  when  they  state  “the  approximation  at  resolution  J  depends  only 
on  the  scaling  function,  and  not  on  the  corresponding  wavelets”[18]. 

4.3  Entropy-Based  Best  Basis  Selection 

Coifman  and  Wickerhauser  [11]  developed  a  technique  for  finding  the  “best”  basis  for  a  given  signal 
using  orthonormal  wavelets.  They  create  a  library  of  known  orthonormal  wavelet  bases,  B.  The  best 
basis  is  that  family  of  wavelets  from  B  that  minimizes  the  entropy  of  the  discrete  decomposition.  Given 
a  vector,  fk,  the  A^-level  discrete  decomposition  consists  of  N  detail  vectors,  for  j  =  1, 2, ...  TV. 
The  entropy  of  the  decomposition  is  defined  as  [1 1] 

«  =  (422) 

3 

However,  because  the  search  of  all  possible  orthonormal  bases  in  a  complete  library  is  too  computation¬ 
ally  intensive,  Coifman  and  Wickerhauser  restrict  the  library  to  only  rapidly  computable  orthonormal 
bases,  like  the  local  trigonometric  basis  and  the  Haar  basis[ll].  The  basis  calculated  is  only  “best”  in 
the  context  of  their  previously  designed  library. 
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4.4  Matching  Pursuit  with  Time-Frequency  Dictionaries 


Mallat  and  Zheng  [25]  and  Chen  and  Donoho  [9]  identified  a  need  to  generalize  the  orthonormal  wavelet 
decomposition  work  of  Mallat  [23].  Mallat  and  Zheng  point  out  that  a  single  wavelet  basis  is  not  flex¬ 
ible  enough  to  represent  a  complicated  non-stationary  signal,  primarily  those  with  narrow  band  high- 
frequency  support  [25].  They  say  that  it  is  much  like  a  person  trying  to  communicate  with  a  very  lim¬ 
ited  vocabulary.  When  the  need  for  a  word  arises  that  is  not  in  the  vocabulary,  many  other  words  are 
needed  to  compensate  for  the  lack  of  the  right  word.  They  go  on  to  say  that  it  is  difficult  to  detect  a  pat¬ 
tern  from  its  projection  coefficients  because  the  information  is  spread  out  over  the  whole  basis.  Each  of 
these  limitations  will  be  specifically  addressed  in  the  design  approach  of  Chapter  5. 

In  order  to  increase  the  flexibility  in  signal  representation,  Mallat,  et.  al.,  construct  a  dictionary  of 
time-frequency  atoms  of  the  form. 


S,(x)  = 


(4.23) 


where  ||p||  =  1,  f  g(x)dx  ^  0  and  p(0)  ^  0.  Notice  that  (4.23)  looks  very  much  like  the  wavelet  ex¬ 
pression  given  by  Morlet  and  Grossman,  except  that  it  includes  a  phase  modulation  term.  This  term  can 
be  used  to  translate  the  frequency  spectrum  of  g~f{x)  up  and  down  the  frequency  axis,  thereby  allowing 
a  narrowband  spectrum  to  reside  at  high  frequencies. 

Mallat  and  Zheng  show  that  the  dictionary  consisting  of  the  time-frequency  atoms  of  (4.23)  is  com¬ 


plete,  that  is,  that  the  closed  linear  span  of  the  dictionary  vectors  is  ^^(3?).  The  objective  is  to  find  the 
best  representation  of  a  signal  using  the  vectors  in  the  dictionary.  Given  a  signal,  /  e  L^(3?),  its  first 
order  representation  is  given  by 


f  —  if i9'fo)9'fo 


(4.24) 


where  5^0  is  the  time  frequency  atom  such  that  |(/,it7o)P  is  a  maximum.  The  first  term  of  (4.24)  repre¬ 
sents  the  orthogonal  projection  of  /  onto  the  space  spanned  by  g-y^.  Therefore,  (/,  570)570  Rf  and 


Wff  =  \\{f,9jo)9yof  +  \\Rff 
=  \{f,9-ro)f  +  \\Rff 


(4.25) 
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since  ||  =  1.  Now  the  process  can  be  repeated  for  Rf, 

Rf  =  {Rf,  )57i  +  -R  V  (4.26) 

where  the  energy  decomposition  of  /  is  now 

ll/f  =  \{f,9'to)f  +  l(/r37i)P  +  (4.27) 

Each  decomposition  of  the  residual  is  done  by  finding  the  time-frequency  atom,  g-y  that  maximizes  the 
energy  of  the  projection  coefficients,  thereby  assuring  rapid  conveigence.  In  applications  like  pattern 
recognition,  signal  reconstruction  and  enhancement,  the  atoms  chosen  and  their  corresponding  projec¬ 
tion  coefficients  can  be  used  to  characterize  or  represent  the  signal  /. 

However,  the  representation  is  only  as  good  as  the  collection  of  atoms  in  its  dictionary.  Each  g-y  must 
either  be  designed  or  discovered  before  ever  proceeding  with  the  adaptive  algorithm.  Furthermore,  the 
algorithm  will  do  best  if  the  atoms  in  the  dictionary  look  like  the  signal  under  analysis  or  its  components. 

4.5  Multiresolution  Analysis-type  Wavelets 

Abry  and  Aldroubi  [2]  developed  an  algorithm  for  designing  wavelets  for  semi-orthogonal  multiresolu¬ 
tion  analyses.  In  a  semi-orthongonal  MRA,  the  wavelet  vector  spaces,  Wj,  are  orthogonal  to  one  another, 
but  'tpj^k  does  not  form  an  orthonormal  basis  of  Wj,  that  is, 

‘^e,m)  =  ^j,t  •  Pk,m  (4.28) 

where  -  k)  and  /  Sk,m-  Just  as  a  non-orthogonal  scaling  function  can  be 

orthogonalized  using  (3.42),  Abry  and  Aldroubi  showed  that  an  orthogonal  scaling  function  can  be  de- 
orthogonalized  by  an  admissible  discrete  sequence,  such  that  the  result  generates  the  same  MRA  as  the 
original  scaling  function. 

Let  <t)±{x)  be  an  orthonormal  scaling  function  that  generates  the  MRA,  {Vj}.  Let 

00 

o,k(f>j.{x  -  k)  (4.29) 

k=—oo 
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be  a  linear  combination  of  the  orthonormal  basis  of  Vq.  Then,  ^{x)  generates  the  same  MRA,  {Vj}, 
which  can  be  easily  shown  in  the  frequency  domain  [2].  The  Fourier  transform  of  (4.29)  gives 


$(u;)  =  2l(a))$i(a;). 


(4.30) 


Notice  that  from  (3.42),  |A(a;)|  must  by  necessity  be  the  square  root  of  the  Poisson  summation  of  $(a;), 
that  is. 


I^Ml  =  1  |^(a;H-27rfc)p  I  . 

<k—~oo 


(4.31) 


Since  (a:)  generates  an  MRA,  it  must  satisfy  its  2-scale  relation,  (3.88),  which  in  the  frequency  domain 

is  given  by 

(4.32) 

Substituting  (4.30)  into  (4.32)  gives 


'  '  A(^)  1.2 


=  H 


a;\  -  /a;\ 

2/  v^y 


(4.33) 


where  H{u})  is  27r  periodic  and  0  <  A  <  A{uj)  <  5  <  oo  is  the  admissibility  condition  on  the 
sequence,  a*  [2].  Since  ^{x)  satisfies  its  2-scale  relation,  it  too  generates  the  same  MRA  as  ^xl^:)-  The 
new  2-scale  relation  is  given  as 

00 

hk$i2x  -  k)  (4.34) 

k=—oo 

where  hf.  is  the  inverse  Fourier  transform  of  H{ix)  in  (4.33). 

The  same  approach  can  be  applied  to  the  wavelet  basis  of  Wj.  Let  'ip_i{x  -  k)  be  the  orthonormal 
basis  of  Wq  and  let  a  new  wavelet  be  constructed  in  Wq  by  a  linear  combination  of  'il)±_{x  —  k),  then. 


-  k) 

k—~oo 


(4.35) 


given  in  the  frequency  domain  as 


^(w)  =  S(cij)^x(w)- 


(4.36) 
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Since  '0j_(a;)  satisfies  its  2-scale  relation  (3.89)  given  in  the  frequency  domain  by 


’J'i(w)  =  G  (^1)  (I)  ,  (4.37) 

then  substituting  (4.36)  into  (4.37)  gives  [2] 

$(a;)  =  5(a;)G  (1^ 

=  °  (I)  (I) 

where  G(u;)  =  B{u;)G{w/2)  is  the  Fourier  transform  of  the  new  generating  sequence,  gk.  The  new 
2-scale  relation  for  -^(a:)  is  therefore  given  as 

OO 

=  IZ  5/fc0±(2a:  -  k)  (4.39) 

k=—oo 

where  if){x)is  expressed  as  a  linear  combination  of  the  orthonormal  basis  of  Vi .  However,  -ipix)  can  also 
be  expressed  as  a  linear  combination  of  any  basis  of  Vi.  Substituting  for  $j.(u;)  in  (4.38)  using  (4.30) 
gives 


4^(0;) 


The  new  2-scale  relation  for  ‘ip{x)  in  terms  of  ^(x)  is  given  as 

00 

=  IZ  -  k) 

fc=-00 


(4.40) 


(4.41) 


where  gl  is  the  inverse  Fourier  transform  of  (w)  =  B{2oj)G{uj) /A{u)  [2]. 

While  ^  and  i>  are  not  orthonormal,  the  new  generating  sequences,  h  and  g,  can  still  be  used  for 
synthesis  in  the  DWT  algorithm  of  Chapter  3.  However,  the  analysis  filters,  h  and  g,  must  be  found 
from  the  duals  of  ^(x)  and  ^(x),  respectively.  The  Fourier  transforms  of  the  duals,  ^  and  ■0  are  given 
as 


$(a;)  = 


j»(u;) _ 

Zr=-oom^  +  2nk)\^ 


1^(0;)  = 


^(a>) 


(4.42) 

(4.43) 
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so  that 


{^j,ki  ^j,m)  —  ^k,m  (4.44) 

{'^j,ki'^j,m)  —  6k, m-  (4.45) 

Notice  that  (4.42)  and  (4.43)  have  the  same  form  as  (4.30)  and  (4.36)  where 


Er=-ool^(^  +  27r/c)P 

(4.46) 

B(u))  -  ^ 

Er=-ool’*'(w  +  27rA:)p- 

(4.47) 

The  analysis  filters,  h  and  g,  can  now  be  found  using  the  same  technique  as  described  above. 

In  order  to  use  the  relationships  given  above  to  design  matching  wavelets,  one  must  be  able  to  find 
the  appropriate  discrete  sequence,  b^.  Abry  and  Aldroubi  find  by  projecting  the  desired  signal,  f{x) 
into  Wo.  The  resulting  approximation  is  given  as 

00 

/(®)  =  £  h'^{x-k)  (4.48) 

k~—oo 

where  the  discrete  sequence,  bk  are  simply  the  projection  coefficients.  One  drawback  to  this  algorithm 
for  finding  matched  wavelets,  is  that  the  MRA  and  therefore  (/)j_  and  'ipj_  must  already  exist.  Furthermore, 
the  quality  of  the  match  will  depend  on  how  much  the  desired  signal  “looks”  like  the  wavelet  basis  in 
the  first  place.  There  is  still  a  need  to  match  an  orthonormal  wavelet  directly  to  a  desired  signal,  with  no 
presuppositions  on  its  shape. 


4.6  Biorthogonal  Wavelets:  The  Lifting  Scheme 

Sweldens  [32]  developed  a  scheme  for  characterizing  all  biorthogonal  filter  sets  that  are  derivable  from 
any  existing  biorthogonal  filter  set.  A  biorthogonal  scaling  function,  ^{x),  generates  an  MRA  {Vj}  and 
has  as  its  associated  wavelet,  ip{x).  The  dual  scaling  function,  ^(x)  generates  a  different  MRA,  {Vj}, 
and  has  as  its  associated  wavelet,  ^(x)  such  that, 

=  0  (4.49) 
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ihk,  (l>j,e)  =  i’j,e)  =  h,i-  (4.50) 

Each  of  the  functions,  <f),  if;,  and  satisfy  their  own  2-scale  relation  of  the  form  (3.88)  and  (3.89) 
where  the  respective  generating  sequences,  h,  g,  h,  and  g,  are  related  by  the  following  expressions. 


G{u)  =  +  tt) 

(4.51) 

G{uj)  =  e-^H{u  +  tt) 

(4.52) 

where  G{u}),  H{lo),  and  G{(j)  are  the  Fourier  transforms  of  h/.,  gk,  K,  and  gk,  respectively. 

Furthermore,  in  order  to  guarantee  perfect  reconstruction, 

H{ui)H{uj)  +  H{oj  +  n)H{uj  +  n)  =  1.  (4.53) 

Sweldens  proves  that  for  a  fixed  (f>  and  its  associated  h,  if  there  are  two  finite  dual  filters,  h  and  h°, 
both  biorthogonal  to  h,  then  the  two  dual  filters  are  related  by  the  following  frequency  domain  expres¬ 
sion: 

Hiuj)  =  H°{u)  +  -b  7r)5(2cu)  (4.54) 

where  ^(a;)  is  a  trigonometric  polynomial  [32].  Sweldens  goes  on  to  prove  the  converse,  that  if  one  of 
the  dual  filters  is  biorthogonal  to  h  and  they  are  related  by  (4.54),  then  the  other  is  biorthogonal  to  h  as 
well. 

Therefore,  (4.54)  can  be  used  to  design  new  biorthogonal  filters  from  existing  ones  where  the  degrees 
of  freedom  in  the  design  reside  in  the  trigonometric  polynomial,  S{uj).  The  design  technique  given  in  the 
following  theorem  is  called  the  ‘Lifting  Scheme”  because  it  lifts  the  properties  of  one  pair  of  biorthog¬ 
onal  filters  in  order  to  find  another  [32]. 

Theorem  3  (Lifting  Scheme)  [  32]  Take  an  initial  set  of  biorthogonal  filters  {h,  hP,g°,  g}.  Then  a  new 


set  of  biorthogonal  filters  {h,  h,g,g}  can  be  found  by 

H{lc)  =  -f-  G(cj)5'(2tu) 

(4.55) 

G(w)  =  G°(a;)  -  HiLj)S{2w). 

(4.56) 
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The  proof  of  this  theorem  can  be  checked  quite  easily.  Since  {/i,  ^]-  are  a  biorthogonal  set  of 

filters,  they  must  satisfy  (4.51),  (4.52)  and  (4.53).  Substituting  for  and  in  (4.52)  and  (4.53)  using 
(4.55)  and  (4.56)  gives  expressions  for  {h,  h,g,g}  that  satisfy  (4.51)-(4.53). 

Theorem  3  can  be  applied  to  a  wavelet  by  substituting  (4.55)  and  (4.56)  into  the  2-scale  relations  for 
the  wavelet  in  the  frequency  domain. 


’5'(a;)  =  G 


if.  I  ^ 


DKI) 


(4.57) 


which  in  the  space  domain  is  given  as 


00  c>0 

^(®)  =  2  E  Pfc0(2rc-A:)-  E  Sk4>{x-k).  (4.58) 

k-~oo  k=~oo 

This  expression  for  'i(^{x)  can  be  used  to  find  the  optimal  with  respect  to  the  error  between  '^(x)  and 
a  desired  signal,  / (a;).  Once  again,  however,  one  must  have  a  pre-defined  biorthogonal  filter  set  before 
proceeding  with  the  adaptive  design.  Furthermore,  for  a  set  of  biorthogonal  filters,  the  family  of  filters 
that  can  be  “lifted”  is  not  complete,  so  the  choice  of  biorthogonal  filters  from  which  one  starts  is  very 
important. 
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Chapter  5 


Matching  a  Wavelet  to  a  Signal 


Using  the  mathematical  construct  of  the  orthonormal  multiresolution  analysis  (MRA),  an  algorithm  for 
finding  a  wavelet  matched  to  a  desired  signal  will  be  developed.  The  technique  performs  the  matching 
in  the  frequency  domain  on  the  magnitude  and  phase  independently  of  one  another.  This  technique  is 
different  from  the  techniques  outlined  in  Chapter  4  because  it  presumes  nothing  about  an  existing  MRA, 
or  existing  wavelets.  Instead,  it  assumes  the  wavelets  to  be  bandlimited,  and  uses  the  conditions  on  an 
orthonormal  MRA  and  the  desired  signal  itself  to  find  the  “closest”  orthonormal  wavelet.  However,  once 
the  algorithm  is  developed  for  the  orthonormal  MRA  case,  it  will  be  shown  that  the  bandlimit  conditions 
can  be  relaxed  so  that  the  matched  wavelet  no  longer  generates  an  orthonormal  MRA,  but  rather  a  frame, 
thereby  providing  much  more  flexibility  in  its  application. 


5.1  Motivation:  Signal  Detection 

Wavelet  transforms  applied  to  multiresolution  analyses  of  signals  produce  outputs  similar  in  theory  to 
those  of  matched  filters[33].  In  order  to  maximize  the  output  of  the  matched  filter  bank,  it  is  necessaiy 
to  design  a  wavelet  that  “matches”  the  signal  of  interest.  This  approach  can  be  justified  by  applying 
matched  filter  theoiy  to  the  wavelet  decomposition  equations.  The  projection  equation  for  the  detail 
functions,  given  in  (3.30),  is  an  inner  product  integral  and  can  be  rewritten  in  the  frequency  domain  by 
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way  of  Parseval’s  Identity: [10] 


4  =  {f{x),i’j,k)  =  (5.1) 

where  is  the  Fourier  transform  of '4’ j,k{x).  The  energy  of  at  a  par¬ 

ticular  scale,  jo,  and  translation,  ko,  is  given  by  its  squared  magnitude 

l4:i"  =  l(f(‘").»*,t.(2'^))l".  (5.2) 

Applying  Schwarz’  inequality  to  the  right  side  of  (5.2)  gives 

|(F(a;),'^,o,fco(2'“‘^))l''  <  (F(a;),F(a;)){«r^„,fe„(2^V),«',„,fc„(2^'>a;))  (5.3) 

where  the  equality  holds  for 

F(u;)  =  K'J'j„,fc„(2^°u;).  (5.4) 

Therefore,  p  is  maximized  when  f(x)  =  Kip  jo, k^.  Rewriting  (5.4)  in  terms  of  amplitude  and  phase 
gives 

where  dpiio)  and  9-^{uj)  are  the  phase  of  F(w)  and  ^'(a;),  respectively.  If  f{x)  matches  exactly  with 
an  orthonormal  wavelet,  then  =  5k,ko  •  ^j,jo  and  the  decomposition  produces  only  1  coefficient  at 
Oo)  /Jo)[33].  As  the  match  moves  away  from  being  perfect,  the  decomposition  begins  to  distribute  the 
energy  about  (jo,  ^o)  which  according  to  (5.2)  and  (5.3)  is  still  the  maximum.  Given  the  condition  for 
a  maximum  projection  coefficient  at  (jo,  ko)  in  (5.5),  the  problem  at  hand  is  to  develop  a  method  for 
matching  the  complex  spectrum  of  the  wavelet  to  that  of  the  desired  signal  while  maintaining  the  condi¬ 
tions  for  an  orthonormal  MRA.  However,  because  the  conditions  for  orthonormality  are  on  the  spectrum 
amplitude  (Poisson  summation)  only,  the  algorithm  matches  the  spectrum  amplitudes  and  group  delays 
independently.  While  this  approach  is  not  ideal  from  an  optimization  standpoint,  it  will  be  shown  that  it 
still  leads  to  good  matching  wavelets. 

A  difficulty  in  matching  the  wavelet  spectrum  directly  to  that  of  the  desired  signal  is  that  the  condi¬ 
tions  on  an  orthonormal  MRA  are  not  if  and  only  if  conditions.  For  instance,  if  is  orthonormal  and 
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generates  an  MRA,  then  is  also  orthonormal.  The  converse  is  not  true[13][24].  Therefore,  in  devel¬ 
oping  the  matching  algorithm  it  is  necessary  to  propagate  the  conditions  for  an  orthonormal  MRA  from 
the  2-scale  sequence  and  scaling  function  to  the  wavelet,  match  the  wavelet  to  the  desired  signal  under 
those  conditions,  and  then  calculate  the  scaling  function  and  2-scale  sequence  always  guaranteeing  that 
the  conditions  for  an  orthonormal  MRA  are  satisfied. 

5.2  Orthonormal  MRAs 

The  detailed  development  of  orthonormal  wavelet  bases  and  multiresolution  analyses  (MRAs)  in  Chap¬ 
ter  3  will  be  summarized  here  for  convenience.  Recall  that  in  an  orthonormal  MRA,  a  signal,  f{x)  € 
F_i,  is  decomposed  into  a  series  of  detail  functions,  {5j(a;)}  j  =  {0, 1, ,  J},  and  a  residual  low 
resolution  approximation,  fj{x),  such  that 

J 

fix)  =  9jix)  +  fjix).  (5.6) 

j=o 

The  decomposition  is  done  by  projecting  P~^{x)  onto  two  orthogonal  subspaces,  Vj  and  Wj,  where 
and  (-i-)  is  the  direct  sum  operator.  The  projection  produces  p{x)  e  Vj,  a  low  res¬ 
olution  approximation  of  f{x),  and  gj{x)  €  Wj,  the  detail  lost  in  going  from  p-^{x)  to  p{x).  The 
orthonormal  bases  of  Wj  and  Vj  are  given  by  =  2-^l‘^'ip{2~^x  -  k)  and  4)j^k  =  2~il‘^(j){2~^x  -  k), 
respectively;  '(j^ix)  is  the  mother  wavelet  and  ^(a:)  is  the  scaling  function  where 

j  ip{x)dx  =  0  =  0  (5.7) 

J  (j){x)dx  =  1  $(0)  =  1  (5.8) 

and  $(a;)  and  are  the  Fourier  Transform  of  (p{x)  and  respectively.  The  projection  equations 
are  given  as 

OO 

9jix)=  dl2-^{2-^x-k)  (5.9) 

k=~oo 

4  =  (5.10) 

OO 

fjix)=  Y  4‘^~hi‘^~^x-k)  (5.11) 

fc=— OO 
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4  =  {fj-i{x),(f>j,k)  (5.12) 

where  d^.  and  are  the  projection  coefficients  and  (•,  •)  is  the  inner  product.  The  nested  sequence 
of  subspaces,  {T^  },  constitutes  the  multiresolution  analysis. 

In  order  for  the  MRA  to  be  orthonormal:  1)  and  (f)j^k  must  be  orthonormal  bases  of  Wj  and  Vj, 
respectively;  2)  Wj  1  Wk,  for  j  ^  k;  and  3)  Wj  X  Vj,  which  leads  to  the  following  conditions  on  ip 
and  <p  [10][20]: 

=  ^k,m  (5.13) 

{<Pj,k,'>Pj,m)  =0  (5.14) 

~  ^j,£  ■  ^k,m-  (5.15) 


The  Fourier  transform  of  (5.13)  gives  the  Poisson  summation,  which  is  1  for  all  w  when  the  integer  trans¬ 
lates  of  (p{x)  are  orthonormal. 


|$(a; -I- 27rm)p  =  1. 


(5.16) 


Since  (p{x)  €  Fq  C  VLi  and  ip{x)  G  C  F_i,  they  can  be  represented  as  a  linear  combination  of  the 


basis  of  VI 


given  in  the  frequency  domain  by 


(p{x)  =  2  hk<p{2x  -  k) 
k=—oo 
oo 

fPix)  =  2  ^  9k<P{2x  -  k) 


(5.17) 


(5.18) 


(5.19) 


(5.20) 


For  orthonormal  MRAs,  the  sequences  hk  and  Qk  in  (5.17)  and  (5. 1 8)  are  quadrature  mirror  filters  (QMF) 
and  have  the  following  properties  [23]  [35]: 


\H{u)f  +  \G{u)f  =  l  (5.21) 

H{u})H{Lj  +  n)  +  G{u})G{u)  +  TT)  =0  (5.22) 

where  if  (w)  and  G(u;)  are  the  Fourier  transforms  of  hk  and  gk,  respectively,  and  are  therefore  both  27r- 
periodic. 
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5.3  Matching  a  Wavelet  to  a  Signal 


5.3.1  Finding  the  scaling  function  from  a  wavelet 

Before  developing  the  algorithm  for  finding  the  matched  wavelet,  a  means  of  deriving  the  scaling  func¬ 
tion  from  the  mother  wavelet  must  be  derived  [7,  8].  Finding  the  wavelet  from  the  scaling  function  is 
simple  using  (5 . 1 8),  however,  it  is  not  invertible.  In  order  to  derive  an  expression  for  |  $  |  in  terms  of  |  ^  | , 
the  following  conditions  are  required  [20]: 


9k  =  {-ifhi^k 

(5.23) 

$(0)  =  1 

(5.24) 

{4>j,ki  4’j,m)  —  ^k,m 

OO 

(5.25) 

(l>{x)=2  hk<j){2x-k). 

(5.26) 

k=—oo 


Condition  (5.23)  simply  guarantees  (5.22)  to  be  satisfied  always.  Conditions  (5.24)-(5.26)  are  required 
for  ^(x)  to  generate  an  orthonormal  MRA,  thereby  satisfying  equations  (5.13)-(5.22)[20]. 

Substituting  (5.19)  and  (5.20)  into  (5.21)  gives  a  relationship  between  |$(a;)|  and  'J'(c<;)|  in  the  con¬ 
text  of  an  MRA  [13]: 

I^MI^  =  l’®'(2u;)P  -t-  |$(2a;)|^.  (5.27) 

Since  the  matching  algorithm  is  performed  on  sampled  data  (in  the  frequency  domain),  there  is  a  need  to 
develop  an  equation  for  finding  the  sampled  scaling  function  spectrum  from  the  sampled  wavelet  spec¬ 
trum. 


Theorem  4  (Finding  |$(fc)|  from  \^{k)\)  In  an  orthonormal  MRA,  let  $(nAa;)  and  ^(nAw)  be  the 
sampled  scaling  function  and  wavelet  spectra,  respectively,  with  sample  spacing  Au  =  7r/2^.  Any  sam¬ 
ple  0/^1$  I  at  (jj  =  n7r/2^  can  be  expressed  by  the  following  recursive  equation: 


for  n  ^  0 


(5.28) 


which  leads  to  the  following  closed  form  solution 


for  n  ^  0. 


(5.29) 
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Furthermore,  (5.28)  implies 


(5.31) 


|’J'(47rA;)|  =  0  for  all  k  £  Z.  (5.30) 

Proof.  Substituting  cj  =  7rn  in  (5.27)  gives 

^  |$(7rn)p  =  |#(27rn)p  +  |’4'(27rn)p. 

However,  since  X)  |^(w  +  27rn)p  =  1  for  orthonormal  MRAs,  and  $(0)  =  1,  and  ^'(0)  =  0  then 

1  for  n  =  0 

0  for  n  ^  0 

and  (5.31)  can  be  rewritten  as 

f 

1  forn  =  0 

|5'(27rn)|  forn  /  0 

So,  at  integer  multiples  of  tt,  |$|  can  be  computed  directly  from  values  of  Furthermore,  (5.32)  and 
(5.33)  imply  |5'(47rn)|  =  0.  Substituting  for  w  =  ;rn/2  in  (5.27)  gives 


|$(27rn)|  =  < 


(5.32) 


|$(7rn)|  =  .J 


(5.33) 


$ 


/  'im\ 

VTj 


=  |$(7rn)p  +  |'4'(7rn)p  forn  /  0. 


(5.34) 


At  integer  multiples  of  7r/2,  |$|  can  be  computed  from  values  of  |  and  the  previously  calculated  values 

of  |$|.  Repeated  substitutions  leads  to  the  following  closed  form  solution: 


$ 


7rn 


=  E 

p=0 


27rn 

2P 


forn  0. 


(5.35) 


If  |’4'(A:Aa;q,)|  has  a  sample  spacing  of  =  27r/2^,  then  by  (5.29),  |$(A:Aw^)|  has  a  minimum 
sample  spacing  of  Aw^,  =  27r/2^^+i  and  i  can  take  on  values  of 


£={0,1,...,M}.  (5.36) 

□ 

Interestingly,  as  £  — >  oo,  equation  (5.29)  approaches  Daubechies’  result  for  continuous  u  [13], 

OO 

for  a;  /  0.  (5.37) 

p=0 

However,  the  beauty  of  (5.29)  over  (5.37)  is  that  the  discrete  form  yields  an  exact  result  for  a  discrete 
spectrum  with  a  finite  number  of  recursive  steps. 
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5.3.2  Properties  of  the  wavelet  spectrum  amplitude 


The  next  step  is  to  derive  constraints  on  |’J'|  that  are  necessary  and  sufficient  to  guarantee  is  an 
orthonormal  basis  of  T^  [7,  8]. 

Theorem  5  (Guaranteeing  Orthonormality)  The  following  condition  on  |  |  is  both  necessary  and 

sufficient  to  guarantee  that  =  h,m,  where  <pj^k  is  derived  from  |^(a;)|; 

00  00  9 

S  1^  ((2”'‘'^‘^  +  27rm))|  =1.  (5.38) 

m— -00  n— 0 

Proof  See  Appendix  A. 

While  (5.29)  and  (5.38)  provide  a  method  for  deriving  an  orthonormal  function,  <j>{x),  from  a  given 
wavelet,  there  is  no  guarantee  that  (f>{x)  generates  an  MRA[13][24].  The  only  condition  remaining  to  be 
incorporated  is  the  2-scale  relation  in  (5.26)  [8].  In  order  for  (f){x)  to  generate  an  MRA,  it  must  satisfy 
its  2-scale  relation  [12][13][24],  which  given  in  the  frequency  domain  is 

$(w)  =  if  j  .  (5  39) 

Repeated  substitution  of  this  recursive  equation  gives 

=  ,5.40, 

where  H{cv)  is  the  Fourier  transform  of  the  discrete  signal,  hk,  and  is  therefore,  27r-periodic.  Equation 
(5.40)  indicates  that  there  must  be  a  certain  structure  in  (j)  in  order  for  it  to  be  a  scaling  function.  In¬ 
corporating  the  infinite  product  into  the  conditions  developed  thus  far  would  be  very  difficult.  A  simple 
way  to  guarantee  that  (5.40)  is  satisfied  is  to  assume  $(u;)  is  bandlimited  [8].  However,  one  is  not  free 
to  choose  any  bandlimits. 

Theorem  6  (Bandlimited  #)  In  a  multiresolution  analysis,  the  spectrum  of  a  bandlimited  scaling  func¬ 
tion,  $(w),  has  maximum  support  given  by 


(5.41) 
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Proof.  See  Appendix  B. 

From  (5.20),  it  is  clear  that  a  bandlimited  scaling  function  that  generates  an  orthonormal  MRA  gives  rise 
to  a  bandlimited  wavelet. 


Corollary  7  (Bandlimited  '$')  In  an  orthononml  MRA  with  a  bandlimited  scaling  function,  the 
sponding  orthonormal  wavelet  has  a  maximum  bandlimit  of 


corre- 


(5.42) 


Proof  See  Appendix  C. 

It  is  interesting  to  note  that  the  bandlimits  derived  above  are  identical  to  those  of  Meyer’s  wavelet[26]. 
However,  what  has  been  derived  here  are  the  maximum  bandlimits  for  not  only  Meyer’s  wavelet  but  any 
bandlimited,  orthonormal  wavelet. 

In  order  to  complete  the  groundwork  for  the  spectrum  amplitude  matching  algorithm,  a  set  of  equa¬ 
tions  for  sampled  spectra,  similar  to  Theorem  5  and  Corollary  7  is  needed. 


Theorem  8  (Guaranteeing  an  Orthonormal  MRA)  LetY{k)  =  |'If(fcAa;)|^,  k  e  Z,  where  Au  = 
27r/2^.  The  necessary  and  sufficient  condition  on  Y  to  guarantee  that  |$(n)|,  found  in  Theorem  4, 
generates  an  orthonormal  MRA  is  given  as  follows: 


E  E>^ 

m=— oo  p=0 


(n  4-  2^+^m)  )  =  1 


2P 


(5.43) 


iM-l 


/3< 


nM 

—  {n  +  2^+^m) 


<2M+2/3  £  =  {0,1,...,M}. 


Proof  See  Appendix  D. 


(5.44) 


Because  the  wavelets  being  designed  are  assumed  to  be  real,  the  magnitude  of  the  wavelet  spectrum 
is  even,  |f'(a;)|  =  |’®'(-w)|,andonly  the  spectra  for  positive  frequency  indices,  k,  in  the  passband  need 
be  matched.  The  conditions  in  Theorem  8  for  k  >  0  generate  a  set  of  L  linear  equality  constraints  in 
F(fe)  of  the  form 

L 

^Q:ifcF(A:)  =  l  for /c  =  {[2^/3],... ,  [2^^+V3J}  (5.45) 
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where  ai^  G  {0, 1, 2}.  Condition  (5.45)  can  be  expressed  in  vector  notation  as 


AY  =  1 


where  A  is  a  L  x  2^  matrix  given  by 

A  =  {ay  €  {0,l,2};i  =  j 

and  1  is  a  L  X  1  vector  given  by 

i^  =  {i  1  ...  1}. 


(5.46) 


(5.47) 


(5.48) 


5.3.3  Matching  Spectrum  Amplitudes 

By  virtue  of  the  limits  given  in  (5.42),  the  desired  signal  must  be  dilated  in  such  a  way  that  the  energy 
in  this  passband  is  a  maximum.  This  dilated  spectrum,  F{u),  is  the  starting  point  for  the  matching  al¬ 
gorithm. 


Theorem  9  (Matched  Wavelet  Amplitude)  Let  W  and  Y  be  vectors  containing  the  samples  of 
|F(fcAu;)p  and  |'J'(A:Aa;)p,  respectively,  in  the  passband: 


W  =  {|F(fcA6,;)|2;A:  = 

[2^/3],...,  [2^+73]} 

(5.49) 

Y  =  k  = 

[2^/3]  , . . . ,  [2^+73]  } 

(5.50) 

where  F{u)  is  the  spectrum  of  the  dilated  signal  being  matched  and^{u})  is  the  matched  wavelet  spec¬ 
trum.  If  the  error  to  be  minimized  is  given  by 


(Iw-Yr(lw-Y) 

then  the  optimal  wavelet  power  spectrum  is  given  by  the  following  expression 


(5.51) 


Y  =  -W  -f  A^(AA^)-i(l  -  -AW) 

Cl  a 


(5.52) 


where 

_  l^(AA^)-iAW 
““  l^(AAr)-ii 


(5.53) 
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and  (AA^)  is  full  rank.  The  match  error  is  given  by 


(1  -  ^AW)^(AA^)-^(1  -  ^AW) 


4w^w 


(5.54) 


The  resultant  wavelet  is  orthonormal,  and  the  scaling  function  it  generates  by  way  of  (5.29)  generates 
an  orthonormal  MRA. 


Proof.  See  Appendix  E. 

Notice  that  the  error,  E,  in  (5.54)  has  the  form  of  a  Mahalonobis  distance,  where  (W^W /a^)AA^  acts 
like  a  covariance  matrix  [28].  This  implies  that  the  solution  is  “closest”  to  the  desired  signal  spectrum 
where  the  distance  measure  is  given  in  (5.51).  The  error  in  the  match  is  a  function  of  the  deviation,  or 
distance,  of  AW  from  1.  If  AW  is  already  1,  then  the  Poisson  summation  of  F{u})  is  1,  f{x)  is  an 
orthonormal  wavelet  and  (f){x),  calculated  in  (5.29),  generates  an  orthonormal  MRA.  As  F{u)  moves 
away  from  being  orthonormal,  AW  moves  away  from  1  and  the  error  of  the  match  increases.  However, 
the  resultant  wavelet  is  still  the  closest  [minimum  distance  as  defined  by  (5.51)]  orthonormal  wavelet  to 

fix)- 

Note:  This  amplitude  matching  algorithm  produces  a  sampled  wavelet  that  in  the  discrete  domain 
satisfies  all  the  properties  of  an  orthonormal  wavelet  While  the  Poisson  summation  for  the  resultant 
scaling  function  and  wavelet  are  shown  to  be  1  at  the  samples  ofu,  there  can  be  nothing  said  about  the 
values  of  the  Poisson  summation  between  the  discrete  samples.  Therefore,  the  resultant  wavelet  spec¬ 
trum  from  Theorem  9  is  a  discrete  approximation  to  a  continuous  orthonormal  wavelet 

The  first  half  of  the  problem  posed  by  (5.5)  has  been  solved,  that  of  finding  the  optimal  wavelet 
spectrum  amplitude  with  respect  to  the  input  spectrum.  The  next  two  sections  develop  the  algorithm  for 
matching  the  wavelet  phase  to  that  of  the  desired  signal. 


5.3.4  Properties  of  the  wavelet  spectrum  phase 

It  would  be  convenient  to  simply  set  the  phase  of  to  the  phase  of  the  desired  signal  spectrum,  F, 
thereby  cancelling  the  complex  exponentials  in  (5.5).  However,  just  as  in  the  previous  section  where 
it  was  shown  that  has  specific  constraints  on  its  amplitude,  here  it  will  be  shown  that  ^  has  specific 
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constraints  on  the  structure  of  its  phase  as  well  [8],  By  repeated  substitutions  of  the  equations  in  (5.19) 

and  (5.20)  and  the  Founer  Transform  of  (5.23),  G{cv)  =  +  tt),  one  gets  the  following  infinite 

products  [10][12][23], 

where  H{u))  is  27r-periodic.  Taking  the  phase  of  both  sides  gives 

(5.55) 

(5.56) 

CX)  ,  . 

m=l  / 

(5.57) 

(5.58) 

where  0$,  6^,  and  dp  are  the  phases  of  and  H,  respectively,  and  6h{(jo)  is  27r-periodic.  The  deriva¬ 

tives  of  (5.57)  and  (5.58)  give  the  negative  of  their  group  delays. 

00  . 

(5.59) 

(5.60) 

where  A$(a;)  d0$(a;)/dw  and  —  d^5-(c<;)/da;,  and  A(a;)  =  d0/r(w)/c?w  is  27r-periodic.  Equation 

(5.60)  can  be  simplified  further  by  letting  ra,(u;)  =  Avp  +  1/2  so  that 


The  wavelet  s  group  delay,  r$(w),  will  be  matched  to  the  group  delay  of  the  desired  signal,  rir(6<;). 
Before  proceeding,  it  is  important  to  note  some  of  the  properties  of  the  group  delays  of  $,  and  H. 

Theorem  10  (Properties  of  A$,A^  and  A)  Lei  A$(u;)  =  (i0$(u;)/da;  W  A^M  =  de^{u;)lduj, 

where  eq,{uj)andd^{u)  are  the  phase  functions  of  ^u})and^{u),  respectively.  LetX{u))  =  deH{uj)/duj 

where  is  the  27r-periodic  phase  ofH{u;).  Then  A^{u;),  A^(a;)  andXiu)  have  the  following  prop- 
erties: 

A$(a))  =  A$(-a;)  (5  62) 
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(5.63) 


A^(a;)  =  Aif  (— w) 


and  therefore. 


A(a;)  =  A(— w) 

\  _  _ 

K^{Lo)du  =  —  X{ijo)du;  =  A  where  X  E  ZS 

-00  -^TT  y_7r 

fOO  I 

/  A^;(a;)da;  = 

J~oo  Z 


/oo 

=  0. 

-00 


(5.64) 

(5.65) 

(5.66) 


(5.67) 


Proof.  See  Appendix  F. 

The  mean  group  delay  of  a  signal  provides  information  about  the  translation  properities  of  that  signal. 
For  instance,  the  mean  group  delay  of  H,  X,  can  take  on  values  of  A  =  ±m  for  me  Z,  since  hk  can  only 
shift  in  unit  increments.  Likewise,  ^(x)  can  only  shift  by  unit  increments,  since  its  mean  group  delay  is 
also  A.  However,  •0(x)  is  fixed  on  the  x-axis  since  its  mean  group  delay  is  a  constant  and  independent 
of  A,  that  is,  shifts  in  hf..  These  same  shifting  relationships  between  </>(x)  and  can  be  derived  from 
(5.17)  and  (5.18). 


Matching  the  group  delay  of  a  desired  signal  to  the  group  delay  of  a  wavelet  given  in  (5.61)  cannot 
be  done  in  the  same  manner  as  the  amplitude  matching  since  there  are  additional  periodicity  constraints 
on  A(a;).  Furthermore,  there  is  still  the  problem  of  finding  the  phase  of  $  from  the  phase  of  To  solve 
both  problems,  one  period  of  A(a;)  is  modeled  as  a  polynomial  of  order  R  [8].  Because  A(a;)  is  an  even 
function  (Theorem  10),  the  polynomial  has  only  even  exponents 

R/2 

=  £cva;2''n(|^)  (5.68) 

r=:U 


where  Cr  are  the  coefficients 


of  the  polynomial  and  n(a;)  is  the  “rect”  function  defined  as 


n(a;) 


1  2  —  ^  ^  2 
0  elsewhere. 


Now  express  A{c<;)  as 


(5.69) 


OO 

\{u,)  =  Xxiiv  —  2'Kk) 

fc  — —  00 
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(5.70) 


oo  till 

k=-oor=0 

Like  the  amplitude  matching  algorithm,  the  phase  matching  algorithm  is  developed  for  discrete  samples 

of  the  spectrum.  Let  Ao;  =  27r/T,  and  N  be  the  number  of  samples  in  F{nAu}).  Equation  (5.70)  can 

/ 

be  rewritten  in  discrete  form  as 

=  E  ^)  (5.71) 

r=0  k=-P/2  ^ 

where  P  =  N/T  is  the  number  of  periods  over  N  points  and  -iV/2  <n  <  N/2.  The  discrete  form  for 
A(n)  can  now  be  written  in  vector  notation  as 


A  =  Be  (5.72) 

where  A  is  a  iV  x  1  vector,  B  is  a  iV  x  {R/2  +  1)  matrix  and  c  is  a  {R/2  +  1)  x  1  vector.  The  elements 
of  B  are  given  as 

77  —  kT 

E  (5.73) 

k=-P/2 

Substituting  (5.72)  into  (5.59)-(5.61)  gives  a  matrix  equation  for  A$,  A^, 


and 


A<^  =  D$c 

(5.74) 

Acj  ^ 

=  — ^  +  D^,c 

(5.75) 

=  D^c 

(5.76) 

where 

00 

=  E  2-™B^  (5.77) 

m—1 

and 

^  00 

Dq,  =  ^  2-”^B^  (5.78) 

m—2 

and  the  elements  of  B  q+r  andB_L  are  given  in  (5.73)  wheren  =  (g  +  T)/2andra  =  g/2'”,  respec- 
tively. 
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5.3.5  Matching  Spectrum  Phase 


In  this  section,  the  expression  for  will  be  derived  such  that  it  is  closest  to  the  desired  signal’s  group 
delay.  Fir,  in  a  least  squares  sense  and  yet  has  the  structure  of  (5.76)  and  (5.78).  Let  7  be  the  unweighted 
error  to  be  minimized. 


N/2-\ 

1=  Y.  (rF(n)-r^(n))2. 


(5.79) 


n=-N/2 


Since  the  wavelet  phase  need  only  match  that  of  the  desired  signal  in  the  passband,  weight  the  error 
function  by  a  normalized  weighting  function.  Let  Q{n)  =  Y  (n)  /^Y  (n),  where  Y  (n)  are  the  elements 

of  Y  derived  in  Theorem  9.  The  weighted  error  function  becomes 

N/2-1 

7n  =  E  mn){TFin)-U{n))]'^ 

n=-N/2 


N/2-1 

=  E 

n=-N/2  I 


R/2 


t2 


Cl{n){T  p{n)  E!  (^ydn.r) 


r=0 


(5.80) 


where  {dn,r}  are  the  elements  of  D4,  in  (5.78).  Rewriting  (5.80)  in  vector  notation  gives 


7  =  (f  F  -  D,iyc)’’(FF  -  D^c) 


(5.81) 


where  the  elements  of  Fp  are  the  non-zero  values  of  {fi(n)rF(n)}  and  the  elements  of  are  the 
corresponding  non-zero  values  of  {n(n)d„,r}-  The  vector,  c,  which  minimizes  7  is  found  by  setting 
Vc7  =  0,  giving 

c  =  (D|D,p)-iD^fF  (5.82) 

where  is  full  rank.  It  follows  that  the  group  delay  of  the  wavelet  can  be  found  by  substituting 

(5.82)  into  (5.76), 


F —  IDi^c. 


(5.83) 


Since  c  in  (5.82)  gives  the  best  estimate  of  c,  it  can  be  substituted  into  (5.72),  (5.75)  and  (5.74)  to  cal¬ 
culate  the  best  estimates  of  A,  and  A^,, 


A  =  Bc 


A(j 


A$  =  ^Dq(C  —  —  2 


(5.84) 
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(5.85) 


A$  =  D$c  —  D$c 


(5.86) 


where  D«,c  and  D$c  are  the  means  of  and  D$c,  respectively.  The  means  in  (5.85)  and  (5.86)  are 
subtracted  so  that  A;|r  and  A$  have  the  properties  of  Theorem  10.  Both  Axj,  and  A$  can  be  summed  to 
obtain  the  discrete  phases  of  #  and  'J'  that  when  combined  with  the  magnitudes  from  Theorem  9  give 
the  full  estimate  of  $(nAa;)  and  ^{nAu)  which  satisfy  all  conditions  for  an  orthonormal  MRA.  The 
QMF  filters,  h  and  g,  corresponding  to  the  matched  wavelet  and  its  scaling  function  can  be  found  using 
(5.19)  and  (5.20)  and  the  inverse  Fourier  Transform.  A  flow  chart  of  the  complete  algorithm,  including 
both  the  amplitude  and  phase  matching  algorithms,  is  given  in  Figure  5. 1 . 


5.4  Generating  Matched  Wavelet  Frames 


The  wavelet  design  algorithm  of  Section  5.3  is  based  on  the  conditions  for  an  orthonormal  MRA  and 
as  a  result,  finds  wavelets  such  that  their  associated  scaling  functions  generate  an  orthonormal  MRA. 
The  condition  in  Section  5.3  that  guarantees  an  MRA  is  the  bandlimit  condition  on  $(u;)  given  in  Theo¬ 
rem  6,  which  leads  to  the  the  bandlimit  conditions  on  ^'(aj)  of  Corollary  7.  It  is  important  to  note  that  the 
matching  algorithm  given  in  Sections  5. 3. 3-5.3. 5  does  not  depend  on  the  specific  bandlimits,  but  only 
the  fact  that  they  exist  (that  ‘ip{x)  is  bandlimited). 

This  section  shows  that  if  the  bandlimits  of  the  wavelet  spectrum  are  expanded  beyond  that  of  Corol¬ 
lary  7  and  if  6  €  5i,  then  the  matching  algorithm  generates  a  dyadic  wavelet.  It  will  also  be  shown  that  if 
b  =  nbo,  that  is,  sampled  with  sample  spacing  6o,  then  the  resultant  wavelet  generates  a  frame  of 
and  there  is  redundancy  in  the  decomposition  of  the  desired  signal  using  the  matched  wavelet. 

Theorem  11  (Matched  Dyadic  Wavelet)  Let  F{u)  be  the  spectrum  of  an  input  signal.  Let  ’5'(u;)  be 
the  spectrum  of  a  bandlimited  wavelet  such  that  ^{u)  0for2Tta  <  \u\  <  2nl3anda  <  1/3, 
(3  >  4/3.  Furthermore,  letW  andY  be  vectors  containing  the  samples  of  \F{kAu})\'^  and  \^{kAu>)\'^, 
respectively,  in  the  passband: 


W=  {|F(A:Au;)|2;fc=  [2^  a],. 

1 - 

to 

(5.87) 

(5.88) 
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Figure  5, 1 :  Algorithm  Flow  Chart 


Then  the  power  spectrum  of  the  optimal  matched  wavelet  is  calculated  using  the  algorithm  in  (5.51)— 
(5.54)  of  Theorem  9.  Furthermore,  the  resultant  wavelet,  given  by, 

where  9p{u))  is  the  phase  of  F{lo),  is  a  dyadic  wavelet  satisfying  the  stability  condition^ 

0  <  A<  Z{()  =  <B<oo  (5.90) 

j=—oo 

and  A  =  B  —  1. 


Proof 

While  ’J'(w)  does  not  satisfy  thebandlimits  of  Theorem  7,  it  does  satisfy  (5.37)  given  in  discrete  form 
in  (5.35),  where  $(w)  is  the  Fourier  Transform  of  some  function,  (j)(x)  (no  longer  a  scaling  function), 
and  as  specified,  $(0)  =  1.  Substituting  a;  =  in  (5.37)  gives 


j=o 

00 


E  |®(n)r. 

j=i-i 


(5.91) 


OO  _ 

=  E  |4'(n)| 

j=~oo 


Taking  the  limit  of  both  sides  as  ^  cx)  gives 

lim  $  (^)  ^ 

i—^OO  \2^/ 

|$(0)|2  = 

1  = 

which  satisfies  the  stability  condition, 

OO 

0<^<  ^  |«'(2^e)|  <B<oo 

j=~oo 


(5.92) 


(5.93) 


where  A  —  B  =  !.□ 

It  is  worth  noting  that  while  the  matched  wavelet  derived  with  Theorem  1 1  does  not  have  an  associated 
scaling  function  and  is  not  associated  with  an  MRA,  it  still  must  satisfy  its  Poisson  summation.  There¬ 
fore,  integer  translates  of  the  matched  dyadic  wavelets  are  orthonormal.  The  violation  of  the  MRA  con¬ 
struct  is  due  to  the  dyadic  wavelets  being  non-orthonormal  with  wavelets  on  a  different  scale,  that  is, 

)  =  Pj,t  ■  4,m  (5 .94) 
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where  pj^i  7^  5j^t. 

The  algorithm  in  this  theorem  is  identical  to  that  of  Theorem  9  with  two  exceptions;  1)  the  bandlim- 
its  on  exceed  those  required  by  Theorem  7;  and  2)  the  phase  of  the  resultant  wavelet  is  simply 
the  phase  of  the  desired  signal.  When  the  bandlimits  are  constrained  to  those  given  in  Theorem  7,  the 
resultant  wavelet  from  the  amplitude  matching  algorithm  is  an  orthonormal  wavelet  whose  associated 
scaling  function  generates  an  orthonormal  MRA.  If  the  bandlimits  exceed  those  of  Theorem  7,  then  the 
matched  wavelet  satisfies  the  stability  condition  (5.90)  and  is  therefore  a  dyadic  wavelet.  Because  the 
resultant  wavelet  has  no  associated  scaling  function,  it  does  not  have  to  exhibit  the  phase  structure  re¬ 
quired  of  wavelets  in  an  orthonormal  MRA  and  can  therefore,  be  any  function.  In  order  to  maximize  the 
matched  filter  output,  the  phase  of  the  wavelet  is  set  to  that  of  the  desired  signal. 

Decomposition/reconstruction  of  a  signal  using  a  dyadic  wavelet  requires  the  shift  paramter,  b,  to 
be  continuous.  In  practice,  b  is  discrete.  Kaiser  [19]  (Theorem  6.1)  proves  that  bandlimited  wavelets, 
critically  sampled  in  b  generate  frames  of  L^(3?). 

Theorem  12  (Bandlimited  Wavelets  and  Frames)  Letip{x)  be  a  real  valued,  bandlimited  wavelet  and 
be  its  Fourier  Transform.  Then  the  family  of  wavelets, 

Ao,j,n  =  2-^{2-jx-nbo)  ,(5.95) 

is  a  frame  of  L?  (3?)  where  60  is  the  sample  rate  parameter  that  satisfies  the  Nyquist  criteria.  The  wavelet 
transform  is  given  by 


Wfi2fn2^bo)  =  if,  Ao,j,n) 

f°°  i 

=  /  f{x)2-^i2-^x-nbo)dx. 

J  —  00 

The  Inverse  Wavelet  Transform  is  given  by 

f  i^)  ~  bo^^^Z  j'ti2^bo)2~2-ip  ^2~^x  —  nbo) 
j  ^ 

where  'ip{x)  is  the  dual  of  tplx)  and  its  Fourier  Transform  is  given  by 


m  = 


(5.96) 


(5.97) 


(5.98) 
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and 


0  <  ^  ^  1$  1^  <  B  <  00.  (5.99) 

3 

Proof 

Let  ^(x)  be  a  real-valued,  bandlimited  wavelet,  then  it  must  satisfy  (3.6H3.9),  Because  has  finite 
support  in  the  right  side  of  (3.9),  where  a  =  2f  can  be  represented  by  its  Fourier  series  expansion, 

2i^{¥0F{0  =  2^bo  (5.100) 

n 

where 

/OO  ■ _ 

-00 

=  Wfi2^,n2jbo)  (5.101) 

which  can  be  shown  by  letting  a  =  2^  and  b  =  n2^bQ  in  (3.8).  Substituting  (5.101)  into  (5.100)  gives 
2t^^(^F(0  =  2^'6o  E  W"/(2^ (5.102) 

n 

Multiplying  both  sides  by  and  summing  over  j  gives 

Y, !«’  (2^'e)  1^  F(0  =  ^'o  E  E  ^/(2^  n2^bo)22  ^  (2^^)  (5  jq3) 

3  j  n 

Again,  lety(^)  ==  Ylj  l'J'(2-'C)P  where  0  <  A  <  y(^)  <  5  <  00,  then 

FiO  =  ^0  E  L  Wf{2^,n2^bo)2h'-^ (2^^) 

j  " 

=  6oEE^/(2^”2J^o)2^^(2^e)e"*^’^^"2''’“.  (5.104) 

3  " 

Taking  the  inverse  Fourier  Transform  of  (5. 1 04)  gives  the  expression  for  the  inverse  Wavelet  Transform 
with  discrete  scale  and  shift. 


where 


□ 


f{x)  =  bo^'Y^Wf{2^,n2^bo)2  ^•ij){2  ^x-nboj  (5.105) 

3  " 


m  = 


E,i’j'(2Jor 


(5.106) 
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Chapter  6 


Examples 


In  this  chapter,  several  examples  are  given  that  demonstrate  the  performance  of  both  the  orthonormal 
and  dyadic  wavelet  matching  algorithms. 

6.1  Orthonormal  Wavelets 

The  performance  of  both  the  magnitude  and  phase  matching  algorithms  will  be  demonstrated  with  sev¬ 
eral  examples.  In  each  of  the  following  examples,  N  —  512,  and  Aoj  =  27r/16  so  that  M  =  4.  In 
each  of  the  figures  shown,  the  input  signal  is  a  dotted  line  and  the  matched  signal  is  a  solid  line.  With 
M  =  4,  the  non-zero  frequency  indices  in  (5.45)  are  fc  =  {6, 7, ... ,  21}.  The  equality  constraints  in 
(5.43)  and  (5.44)  of  Theorem  8  generate  L  =  11  equations  and  16  unknowns  represented  by  the  A  ma¬ 
trix  in  Figure  6.1.  Recall  that  (5.43)  was  derived  by  constructing  |$(A:Aa;)p  from  the  sum  of  repeated 
dilations  |^(fcAa;)pthen  replicating  that  sum  every  27r  and  summing  giving  the  Poisson  summation  for 
#,  which  should  be  1  everywhere.  This  process  is  illustrated  in  Figure  6.2,  taking  into  account  the  ban- 
dlimits  of  (5.44)  and  solving  for  |$(A:Aa;/2)p  instead  of  |$(fcAa;)p.  Dilating  |  t(A:Aa;/2)p  gives  back 
|'I'(A;Aa;)p  which  is  the  first  row  of  Figure  6.2.  Dilating  again  produces  the  second  row  and  so  on.  The 
pattern  exhibited  in  the  first  5  rows  of  Figure  6.2  exist  for  negative  frequency  indices  as  well.  Now, 
replicating  that  structure  by  47r  (since  we  have  just  constructed  |$(A:Aa;/2)p)  is  done  by  shifting  and 
centering  it  about /:  =  32,  since  Aw  =  tt/S,  which  produces  the  structure  to  the  lower  right  of  Figure  6.2. 
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Figure  6. 1 :  Constraint  Matrix  -  A  for  27r/3  <  a;  <  Stt /3 

These  values  were  the  negative  frequency  values  of  the  original  structure.  The  positive  frequency  values 
reside  on  the  opposite  side  of  A:  =  32.  Replicating  any  further  in  either  the  positive  or  negative  direction 
would  not  produce  any  additional  overlap.  Therefore,  the  complete,  unique  set  of  11  independent  equa¬ 
tions  (enclosed  in  the  dotted  box  in  Figure  6.2)  can  now  be  obtained  from  columns  k  =  6, . . . ,  16,  by 
adding  columnwise  and  remembering  that  the  sum  for  any  column  must  be  1.  The  coefficients  of  these 
11  equations,  16  unknowns  are  represented  in  the  condition  matrix  A. 


6.1.1  Meyer’s  Wavelet 


Meyer’s  wavelet  is  a  bandlimited  wavelet  that  forms  an  orthonormal  basis  of  L2(3?).  If  it  is  used  as  the 
desired  signal,  the  matching  algorithm  should  produce  Meyer’s  wavelet  exactly.  Since  the  wavelet  is  real 
and  symmetric,  there  is  no  need  to  match  the  phase.  Meyer’s  wavelet,  shown  in  Figure  6.3,  is  defined 
in  the  frequency  domain  as  [10, 13]: 


< 


0 


M<f  orH>f 

¥  <  kl  <  f 
f  <  kl  <  f 


(6.1) 
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Figure  6.2:  Construction  of  the  constraint  matrix  A 
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where 


0  7  <  0 

t^(7)  =  <  7  0  <  7  <  1  •  (6.2) 

1  7  >  1 

The  matched  wavelet  spectrum  in  the  passband  is  shown  in  Figure  6.4.  The  scaling  parameter,  a,  and 


Figure  6.3:  Meyer’s  wavelet 

the  match  error  were  found  to  be  1  and  0,  respectively,  since  AW  =  1  in  (5.53).  The  scaling  function 
associated  with  Meyer’s  wavelet ,  calculated  using  (5.35)  is  shown  in  Figure  6.5. 
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Figure  6.4:  Amplitude  Match  in  the  passband  -  Meyer 


Figure  6.5:  Meyer’s  scaling  function 
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6.1.2  Gabor’s  Wavelet 


Gabor’s  or  Morlet’s  wavelet  is  also  real  and  symmetric,  but  does  not  generate  an  orthonormal  basis  of 
Gabor’s  wavelet  is  a  modulated  gaussian  given  by  [13]: 

fG{x)=Ke(^')  cos{u;ox)  (6.3) 

and  in  the  frequency  domain  as 

The  parameters,  K,  a,  and  wo  were  chosen  so  that  ||^G(a;)||  =  1  and  ’5'g(w)  was  centered  in  the  pass- 
band  of  the  matching  algorithm.  Figure  6.6  shows  the  Gabor  wavelet  and  Figure  6.7  shows  its  frequency 
spectrum  along  with  its  Poisson  summation.  The  Gabor  wavelet  is  clearly  not  orthonormal.  Figure  6.8 


Figure  6.6:  Gabor’s  wavelet 

shows  the  results  of  the  amplitude  match  in  the  passband.  The  scaling  parameter,  a  =  0.9552  and  the 
match  error  was  0.0179.  The  resultant  wavelet  and  scaling  function  spectra  with  their  Poisson  summa¬ 
tions  are  shown  in  Figures  6.9  and  6. 10,  respectively.  They  are  both  orthonormal.  The  matched  wavelet 
and  the  desired  signal  are  shown  in  Figure  6.1 1  and  the  corresponding  scaling  function  in  Figure  6.12. 
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Figure  6.8;  Amplitude  match  in  the  passband  -  Gabor’s  wavelet 
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Figure  6.11:  Matched  wavelet  vs  Gabor’s  wavelet 


Figure  6.12:  Scaling  function  -  Gabor 
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6.1^  Daubechies’  D4  wavelet 


Daubechies’  D4  wavelet  and  scaling  function  are  shown  in  Figure  6.13.  Let  the  desired  signal  be  the  D4 


\|/(X) 

_ ^ 

Figure  6.13:  Daubechies’  D4  Wavelet  and  Scaling  Function 


wavelet,  given  by  foix).  The  desired  signal  power  spectrum,  {W (n)|n  =  -256, . . . ,  255},  is  given  as 


W{n)  =  { 


\FD{nAoj)\^ 


for|n|  =  {6, 7,  ...,21} 


(6.5) 


10  forlnl  {6,7,...,21} 

where  Fd{uj)  is  the  Fourier  Transform  of  foix).  The  truncated  spectrum  is  shown  in  Figure  6.14  along 
with  its  Poisson  summation.  Notice  that  the  Poisson  summation  is  no  longer  1  everywhere,  due  to  trun¬ 
cation.  Y (k)  is  found  using  (5.52)  and  (5.53)  where  a  =  0.9124  and  W  =  {ir()fc)|A;  =  6, 7, . . . ,  21}. 
The  results  of  the  match  in  the  passband  are  shown  in  Figure  6.15.  The  full  matched  wavelet  spectrum  is 
constructed  by  reflecting  Y (k)  onto  the  negative  axis,  and  taking  its  square  root.  The  matched  wavelet 
spectrum  and  its  Poisson  summation  are  shown  in  Figure  6.16.  The  wavelet  is  clearly  orthonormal.  The 
scaling  function  magnitude  spectrum  is  calculated  using  (5.28)  and  is  shown  in  Figure  6.17  along  with 
its  Poisson  summation.  The  scaling  function  is  orthonormal  as  well.  Because  we  have  chosen  the  ban- 


dlimits  of  Theorem  7,  we  are  guaranteed  that  the  resultant  scaling  function  generates  an  orthonormal 
MRA.  If  we  were  to  stop  here,  we  would  have  a  symmetric  wavelet  that  comes  closest  to  foix)  and  its 
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associated  scaling  function.  The  asymmetry  in  the  D4  wavelet  is  caused  by  the  structure  in  its  phase. 
The  first  step  in  finding  the  matched  phase  is  to  find  the  group  delay  of  the  desired  signal,  Ap,  which  is 
done  using  the  following  process: 

1.  Calculate  Aw)  =  F£)(nAw)/|F£)(nAw)|. 

2.  Interpolate  across  samples  of  F^(nAw)  where  |F£)(«Aw)|  =  0. 

3.  Ap  =  |A^F^(nAa;)|  where  A^  is  the  first  difference  operator. 

Note:  This  procedure  eliminates  the  need  to  unwrap  27r  phase  jumps  caused  by  the  tan~^  operator. 
However,  this  procedure  also  generates  a  group  delay  that  could  have  come  from  one  of  four  phases, 
dp,  —dp,  dp  +  TT,  and  —dp  +  TT.  The  matched  wavelet  is  calcinated  for  each  of  these  possible  phases 
and  then  the  one  closest  to  the  desired  signal  is  chosen. 

Next,  the  matrix  from  (5.78)  is  calculated.  In  these  examples,  N  =  512  and  R  =  16,  mak¬ 
ing  a  512  X  9  matrix.  The  polynomial  coefficient  vector,  c,  is  calculated  using  (5.82)  where 
and  fi?  are  weighted  by  the  normalized  matched  spectrum,  Y,  calculated  above.  Figure  6.18  shows 
Ap  and  A$.  They  match  veiy  closely  since  Ap-  is  the  group  delay  of  a  known  orthonormal  wavelet 


Figure  6.18:  Matched  Wavelet  Group  Delay  vs  desired  -  D4 
and  therefore  must  have  the  proper  structure.  The  group  delay  of  the  scaling  function,  shown  in  Fig- 
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k 

h(k) 

k 

h(k) 

k 

h(k) 

-8 

-0.0082 

-2 

-0.0055 

3 

-0.0407 

-7 

0.0256 

-1 

0.3515 

4 

0.0338 

-6 

-0.0085 

0 

0.5554 

5 

-0.0202 

-5 

-0.0173 

1 

0.2281 

6 

-0.0037 

-4 

0.0298 

2 

-0.0957 

7 

0.0253 

-3 

-0.0322 

Table  6. 1 :  h{k)  for  ip  matched  to  /d 

ure  6.19  with  the  group  delay  of  the  D4  scaling  function,  is  calculated  using  (5.77)  and  (5.86).  The 


Figure  6.19:  Scaling  Function  Group  Delays;  Derived  vs  Truth  -  D4 

matched  wavelet  and  scaling  function  phase  is  found  by  integrating  (or  summing)  A^(n)  and  A$(n). 
The  matched  wavelet  and  its  associated  scaling  function,  shown  in  Figure  6.20,  are  each  found  by  tak¬ 
ing  the  inverse  Fourier  Transform  of  their  complex  spectra.  The  inner  product  of  fo  with  its  matched 
wavelet,  ip,  gives  {. . .  0.0082  -  0.0202  0.9963  0.0292  0.0072  . . .}.  The  resultant  QMF  filters,  h 
and  g,  are  shown  in  Figure  6.2 1  and  the  middle  16  values  of  h  are  given  in  Table  6. 1 . 
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6.1.4  l^ansient  signal 


The  next  example  is  a  transient  sinusoid  given  by  the  following  equation 


frix)  =  a;e“®  cos{2TTfox)u{x)  (6.6) 

where  u{x)  is  the  unit-step  function.  The  transient  signal  in  this  example  was  constructed  by  setting  a  — 
2.0  and  /q  =  0.8,  and  dilating  it  such  that  its  spectrum,  Ft{uj),  had  maximum  energy  in  the  passband, 

27r/3  <  |u;|  <  Stt/S.  Figure  6.22  shows  the  transient  signal  and  Figure  6.23,  its  spectrum  amplitude  and 
Poisson  summation.  The  transient  signal  is  clearly  non-orthonormal.  Following  the  same  sequence  of 

on.! 

0.02 

0.01 

0.00 

4)0! 

4)02 

-0.03 

4)04 

4)05 


A  - 

rf 

Figure  6.22:  Transient  Signal 

steps  described  above  for  the  D4  wavelet,  the  asymmetric  wavelet  that  best  matches  frix)  is  derived. 
Figure  6.24  shows  the  result  of  matching  the  wavelet  in  the  positive  passband,  where  a  =  1.0348.  The 
full  matched  wavelet  and  scaling  function  spectra  along  with  their  Poisson  summations  are  shown  in 
Figures  6.25  and  6.26.  The  group  delays  of  the  desired  signal,  Ap,  and  the  matched  wavelet,  A^, 
are  shown  in  Figure  6.27.  Since  frix)  is  not  a  wavelet,  its  phase  shouldn’t  have  the  required  structure. 
However,  notice  that  the  matched  wavelet  group  delay  does  have  the  required  structure  and  matches  the 
desired  group  delay  very  well  in  the  passband.  The  matched  wavelet  and  its  associated  scaling  function 
are  shown  in  Figures  6.28  and  6.29.  The  inner  product  of  fr  with  its  matched  wavelet,  V'  gives  {.. .  - 
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Figure  6.29:  Scaling  Function  -  transient 

0.0062  -0.0015  0.9834  -0.0412  -0.0935  ...}.  Eventhough /r  is  not  bandlimited,  its  correlation 
with  the  matched  wavelet  still  produces  a  value  very  near  1.0,  with  very  little  spread  in  translation.  The 
resultant  QMF  filters,  h  and  g,  are  shown  in  Figure  6.30  and  the  middle  16  values  of  the  QMF  filter,  h, 
are  given  in  Table  6.2. 
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6.2  Dyadic  Wavelets 


This  section  gives  some  examples  of  optimally  matched  dyadic  wavelets  using  the  algorithm  in  Theo¬ 
rem  11.  For  each  example,  N  =  512  and  Aw  =  27r/16  so  that  M  =  4  as  in  the  previous  section.  The 
bandlimits  will  be'O  <  |w|  <  47r  which  gives  a  constraint  matrix,  A,  that  is  a  16  x  32  matrix  shown  in 
Figure  6.31.  The  frequency  indices  in  (5.45)  are  is:  =  {1, 2, ... ,  32}. 


1  lOIOOOlOOOOOOOlOOOOOOOOOOOOOOI  l" 
OlOIOOOlOOOOOOOlOOOOOOOOOOOOOlOl 
00100100000100000000000100001000 
00010001000000010000000000010001 
00001000010000000001000000100000 
00000100000100000000000101000000 
00000010000001000000000010010000 
00000001000000010000000100000001 
00000000100000000100001000000000 
00000000010000000001010000000000 
00000000001000000000110000000000 
OOOOOOOOOOOIOOOOOOOIOOOIOOOOOOOO 
00000000000010000010000001000000 
0000000000000  1  000  1  000000000  1  0000 
00000000000000101000000000000100 
D0000000000000020000000000000002. 


Figure  6.31:  Constraint  Matrix  -  A  for  0  <  w  <  47r 

6.2.1  Gabor’s  Wavelet 

Gabor’s  wavelet,  fcix),  defined  by  (6.3),  is  dilated  so  that  a  maximum  amount  of  its  energy  in  the 
frequency  domain  is  contained  in  the  expanded  passband  represented  by  the  constraint  matrix,  A  (Fig¬ 
ure  6.3 1 ).  The  result  of  the  amplitude  match  is  shown  in  Figure  6.32.  Notice  that  there  were  samples  that 
were  driven  to  0  by  the  inequality  constraint,  Y(k)  >  0.  As  the  passband  increases,  the  spectrum  of  the 
desired  signal  becomes  closer  to  0,  and  there  is  a  higher  chance  that  the  optimal  solution  would  drive 
the  sample  to  be  negative.  The  inequality  constraint  is  invoked  to  set  that  value  to  0,  thereby  making 
the  solution  sub-optimal.  Because  foix)  is  real  and  symmetric,  there  is  no  phase  to  apply  to  ip{x).  The 
resultant  dyadic  wavelet  is  shown  in  Figure  6.33  along  with  foix).  The  excessive  ringing  at  the  edges 


100 


is  due  to  the  spikes  in  the  spectrum  match  caused  by  the  inequality  constraints.  Figure  6.34  shows  the 
amplitude  spectrum  of  the  matched  dyadic  wavelet  along  with  its  Poisson  summation.  The  algorithm  of 


Figure  6.34:  Matched  Dyadic  wavelet  spectrum  and  poisson  summation  -  Gabor 

Theorem  11,  like  Theorem  9  guarantees  the  Poisson  summation  of  the  resultant  wavelet  to  be  1.  This 
implies  that  integer  translates  of  the  matched  dyadic  wavelet  is  orthonormal! 

6.2.2  Daubechies’  Wavelet 

Daubechies’  wavelet,  shown  in  Figure  6.13,  is  the  last  example  for  the  dyadic  matching  algorithm.  The 
algorithm  should  perform  better  than  in  Section  6. 1  because  the  extended  bandlimits  should  produce  less 
truncation.  Figure  6.35  shows  the  result  of  the  amplitude  match  in  the  passband  and  Figure  6.36  shows 
the  wavelet  spectrum  and  Poisson  summation,  which  is  1  everywhere,  as  expected.  Figure  6.37  shows 
the  resultant  matched  dyadic  wavelet  along  with  Daubechies’  D4  wavelet. 
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Chapter  7 


Applications  to  Image  Object  Detection 


In  this  chapter,  the  matched  wavelet  algorithm  is  used  to  detect  and  identify  an  object  in  an  image.  The 
approach  does  not  use  the  2-D  Discrete  Wavelet  Transform  because  of  its  shift  variance  and  lack  of  ro¬ 
tation  sensitivity  (see  Chapter  3.8).  Instead,  the  2-dimensional  object  is  transformed  into  a  series  of  1- 
dimensional  vectors  using  the  sampled  Radon  Transform  [15].  The  matched  wavelets  for  each  of  the  1  -D 
vectors  is  calculated  and  used  to  detect  the  object,  and  determine  the  aspect  ratio,  scale,  and  angle  of  ro¬ 
tation.  Further  application  of  the  matched  wavelet  algorithm  on  the  object  details  provides  additional 
feature  vectors  that  can  be  used  for  more  detailed  identification. 

This  chapter  is  divided  into  four  sections,  the  first  two  provide  background  on  the  Radon  Transform 
and  image  reconstruction  using  backprojection  and  the  projection  slice  theorem.  Section  7.3  contains 
an  introduction  to  the  Wavelet  Radon  Transform  and  Section  7.4  develops  the  detection  algorithm  and 
gives  results. 

7.1  The  Radon  Transform 

The  Radon  Transform  was  developed  by  its  namesake,  Johann  Radon,  in  his  1917  paper  entitled  “On 
the  determination  of  functions  from  their  integrals  along  certain  manifolds”  [29].  The  Radon  Transform 
allows  one  to  represent  a  2-dimensional  function,  f{x,  y),  by  an  infinite  series  of  1 -dimensional  projec- 
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tions  that  cover  the  x-y  plane.  The  standard  form  of  the  Radon  Transform  is  as  follows: 

qeip)  =  f  fix,y)da  (7.1) 

L 

where  q  is  the  Radon  Transform  of  /,  n  is  the  Radon  Transform  operator,  and  L  is  the  line  defined  by 

p  =  a;cos^  +  j/sin0  (7.2) 

for  p  €  3?  and  0  <  0  <  TT.  As  shown  in  Figure  7.1  the  Radon  Transform  is  formed  by  rotating  the 
coordinate  axis  through  some  continuous  angle,  9,  and  integrating  along  the  new  vertical  axis  (<T-axis). 
The  resultant  transform  is  a  function  of  both  the  new  horizontal  axis,  p,  and  the  angle  of  rotation,  9. 
Using  the  coordinate  transformation  equations  for  a  rotated  coordinate  frame  gives  the  following  form 


Figure  7.1:  Projection  geometry  for  tomographic  processing 
of  the  Radon  Transform, 

qe (p)  =  'R'f  =  j^f{pcos9  -  a  sin 9,psin9  +  a  cos  9)da.  (7.3) 

Figure  7.2  shows  an  image  of  a  2-D  rectangle  function  and  its  Radon  Transform.  Notice  that  at  0°,  the 
Radon  Transform  is  a  “rect”  function  whose  width  is  equal  to  the  width  of  the  rectangular  object.  Like- 
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wise,  at  90°,  the  Radon  Transform  is  a  “rect”  function  whose  width  is  the  same  as  the  height  of  the  rect¬ 
angular  object.  The  Radon  Transform  is  taken  for  0  <  <  tt  because  the  projections  for  tt  <  0  <  27r 


180 

90 

0 

_ ^ 

Figure  7.2:  Radon  Transform  of  a  rectangle  image 

contain  redundant  information.  The  symmetry  relationship  between  projections  taken  from  an  angle 
0  <  0  <  TT  and  an  angle  rotated  an  additional  tt  radians  is  given  as 

qe+Ap)  =  qe{-p)-  (7.4) 

This  relationship  can  be  seen  graphically  in  Figure  7.1.  A  projection  taken  180°  from  the  angle  shown 
will  produce  a  projection  that  is  simply  flipped  along  its  horizontal  axis. 

When  0  is  sampled  with  sample  spacing  A0,  the  resultant  transform  is  referred  to  as  the  Sampled 
Radon  Transform  [15].  The  Radon  Transform  in  Figure  7.2  is  actually  a  Sampled  Radon  Transform 
because  it  was  constructed  with  A0  =  1°. 
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7.2  Image  Reconstruction  and  Backprojection 


Although  the  basic  formulas  for  inversion  of  the  Radon  Transform  were  worked  out  as  early  as  1917 
[29],  very  little  if  any  effort  was  devoted  to  implementing  the  inversion  in  a  practical  situation  prior  to 
the  pioneer  work  in  radio  astronomy  by  Bracewell  [4][15].  The  practical  problem  with  the  inversion  for¬ 
mulas  derived  by  Radon  is  that  they  require  /  (a;,  y)  to  be  a  continuous  function  and  the  Radon  Transform 
to  be  constructed  from  all  lines  in  the  x-y  plane  [15]. 


7.2.1  Reconstruction  by  direct  Fourier  Methods 

Bracewell  [4]  derived  the  relationship  between  the  2-D  Fourier  Transform  of /(x,  y)  and  the  1-D  Fourier 
Transform  of  qeip)  with  respect  to  p.  The  relationship  is  given  as 


(^)  C)  =  Fg{ujx  cos  0  +  ujy  sin 9, 0) 

=  F(7/cos0  — C  sin  sin  Ceos  0) 

/OO 

q9{p)e-^'^^^Pdp  (7.5) 

-OO 

where  rj  and  C  are  the  spatial  frequency  variables  as  shown  in  Figure  7.3,  and  Fg  is  F  rotated  through 
9.  The  result,  later  referred  to  as  the  Projection  Slice  Theorem,  states  that  the  1-D  Fourier  Transform 
of  the  Radon  Transform  for  a  fixed  9  is  equal  to  a  slice  through  the  2-D  Fourier  Transform  of  f{x,  y) 
at  an  angle  9  (Figure  7.3)  [15].  Using  the  Projection  Slice  Theorem,  one  could  approximate  F{u}x,  ojy) 
given  enough  projections  by  forming  a  superposition  of  the  Fourier  Transform  of  each  projection  in  the 
Sampled  Radon  Transform  [5],  given  as 


F{u}x.,UJy)  =I 


(7.6) 


where  Fg^  (y)  is  the  Fourier  Transform  of  the  projection  taken  at  9  =  9i  and  I  is  an  interpolation  oper¬ 
ator  that  resamples  the  data  onto  a  cartesian  grid.  However,  in  practice,  there  is  a  finite  number  of  pro¬ 
jections  and  the  sample  spacing  in  the  frequency  domain  is  nonuniform,  which  causes  spatially  variant 
errors  after  interpolation.  Several  standard  techniques  have  been  developed  to  alleviate  the  non-uniform 
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Figure  7.3:  Projection  Slice  Theorem 

sampling  problem.  Bracewell  and  others,  however,  developed  a  space  domain  solution,  called  backpro- 
jection,  that  works  around  the  non-uniform  sample  grid  in  the  Fourier  domain. 

7.2.2  Reconstruction  by  backprojection 

Instead  of  constructing  F  as  in  (7.6)  and  then  taking  the  inverse  Fourier  Transform,  one  could  take  the 
inverse  Fourier  Transform  of  Fg.  [i],  0)  first,  and  sum  the  results,  given  as 

/(^)2/)  =  ^  qe{xcos9  +  ysin$)d6.  (7.7) 

The  inverse  Fourier  Transform  of  a  slice  in  the  frequency  domain  at  angle,  6,  is  a  2-D  “ridge”  whose  cross 
section  through  the  angle  6  is  fo{p)  and  whose  values  are  constant  perpendicular  to  its  cross  section  [5], 
as  shown  in  Figure  7.4.  The  inverse  Fourier  Transform  of  a  slice  is  referred  to  as  a  backprojection,  since 
it  is  a  1-D  function  smeared,  or  replicated  back  along  its  perpendicular  axis.  The  superposition  of  all 
backprojections  forms  a  layergram  and  gives  a  coarse  approximation  of  the  original  function,  f{x,  y). 
However,  the  impulse  response  of  the  layergram  is  not  a  delta  function.  In  fact,  its  transfer  function 
converges  to  H{y,  0)  =  1/7^  (polar  notation)  as  M  ->  0.  Therefore,  before  taking  the  inverse  Fourier 
Transform  of  the  slice,  it  is  necessaiy  to  multiply  it  by  the  ramp  function,  r{r])  =  \y\,  where  -rj^  < 
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Figure  7.4:  Backprojection  of  qe{p) 

V  <  Vn  and  t]n  is  the  Nyquist  frequency.  The  complete  modified  or  filtered  backprojection  algorithm 
for  image  reconstruction  is  shown  in  Figure  7.5.  The  algorithm  developed  over  the  next  two  sections 
builds  off  of  and  modifies  the  backprojection  algorithm  presented  above. 

7.3  The  Wavelet  Radon  Transform 

Let  Figure  7 .6  represent  a  typical  collection  geometry  for  the  projection  analysis  described  above.  A  2-D 
image,  / {x,  y),  is  represented  by  8  projections,  denoted  as  qi{p),  taken  at  equal  angle  intervals,  where  the 
angle  at  which  the  projection  is  taken  is  di  =  (i  -  1)A0,  where  =  tt/S  and  f  =  1, 2, . . . ,  8.  Taking 
the  continuous  wavelet  transform  (CWT)  of  each  projection  gives  the  Wavelet  Radon  Transform.  Since 
in  practice,  the  input  signals  are  discrete  signals,  the  CWT  is  approximated  using  matrix  multiplication. 
Appendix  G  gives  the  specific  implementation  of  the  CWT  on  sampled  data.  Each  projection  can  now  be 
analyzed  by  way  of  its  wavelet  transform  for  object  detection  and  identification.  For  example,  assume 
a  512  X  512  image  consists  of  a  uniform  background  and  a  rectangular  object  of  width  10  and  height  30 
as  shown  in  Figure  7.7.  Figure  7.8  shows  its  8  projections.  The  projections  have  different  frequency 
content  (some  are  wider,  some  are  narrower),  which  can  be  determined  by  way  of  its  Wavelet  Radon 
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Figure  7.5:  Image  Reconstruction  using  the  Backprojection  Algorithm 


Figure  7.6:  Geometry  for  the  Wavelet  Radon  Transform 
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Figure  7.7:  Image  of  a  rectangular  object 


Transform,  shown  in  Figure  7.9.  Let  the  wavelet  transform  of  qi{p)  be  denoted  as  W^{a,  b),  where  a 
and  b  are  the  scale  and  shift  parameters,  respectively.  In  Figure  7.9,  {a,  b)  is  displayed  at  the  top 

while  Wg  (o,  b)  is  at  the  bottom,  and  for  each  wavelet  transform,  scale  increases  downward.  Notice  that 
the  peak  in  (a,  b)  occurs  at  very  small  scale  and  as  i  increases,  the  peak  moves  to  larger  scales  until 
it  resides  at  the  maximum  scale  at  z  =  5  which  corresponds  to  90°  .  The  trajectory  of  the  peak  contains 
information  about  the  aspect  ratio  of  the  object  and  will  be  used  in  the  next  section  as  a  feature  vector 
for  identifying  objects  in  an  image. 

The  obvious  skew  in  the  CWTs  in  Figure  7.9  is  due  to  the  tl^at  is  required  of  an 

orthonormal  wavelet  (see  Section  5.3.4).  The  phase  term  causes  the  wavelet  to  be  centered  about  x  = 
1/2  and  as  the  wavelet  is  dilated,  the  location  of  its  center  shifts,  which  causes  the  skew  in  the  wavelet 
tranform. 

Notice  also  in  Figure  7.9  that  at  very  small  scales,  the  wavelet  evokes  a  response  from  the  edges  of 
the  rectangle  functions  in  the  projections.  This  is  especially  noticeable  in  projections  3-7.  The  response 
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Figure  7.8:  Projections  of  rectangle  image  at  8  equally  spaced  angular  intervals 


Figure  7.9:  Wavelet  Radon  Transform  of  rectangle  image  using  the  D4  wavelet 
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from  each  edge  (which  is  high  frequency  content  or  small  scale)  moves  toward  the  center  as  the  scale 
increases  until  the  two  intersect  at  the  peak  of  the  wavelet  transform.  At  this  scale,  the  wavelet  is  matched 
to  the  scale  of  the  entire  object,  not  just  the  edges.  Because  the  peak  of  the  wavelet  transform  is  important 
in  analyzing  the  object,  it  is  necessary  that  the  peak  be  maximized,  which  can  be  done  with  the  matched 
wavelet  algorithm  of  Chapter  5. 


7.4  Matched  Wavelets  and  Object  Detection 

The  Wavelet  Radon  Transform  is  a  useful  tool  for  analyzing  objects  in  an  image  where  rotation,  scale, 
and  position  are  unknown.  As  shown  in  Section  5.1,  the  wavelet  transform  is  much  like  a  matched  filter 
and  therefore,  the  response  in  the  Wavelet  Radon  Transform  is  maximized  when  the  wavelets  match  the 
projections  under  analysis.  For  that  reason,  the  wavelet  matching  algorithm  of  Chapter  5  will  be  used  to 
find  the  optimal  wavelets  for  each  projection  of  the  image  object. 

Because  the  optimal  wavelets  are  bandlimited,  they  can  match  only  one  segment  of  a  signal’s  fre¬ 
quency  spectrum  at  any  given  time.  The  matching  algorithm  uses  the  maximum  energy  criterion  for 
choosing  which  segment  to  match.  This  limitation  is  resolved  by  making  multiple  passes  over  the  pro¬ 
jections.  After  each  pass,  a  residual  detail  signal  is  created  by  taking  the  difference  between  the  original 
projection  and  the  detected  signals.  A  second  pass  is  executed  with  the  residual  detail  signals  as  input. 

7.4.1  IVaining  on  a  Known  Object 

Given  an  input  image  containing  a  known  object,  it  is  necessary  to  train  the  detection  algorithm  on 
that  given  object.  Figure  7.10  shows  the  training  procedure.  This  training  procedure  provides  a  set  of 
matched  wavelets  for  two  passes  and  a  maximum  scale  vector  that  provides  aspect  ratio,  location,  and 
rotation  information.  Once  these  functions  and  parameter  values  are  obtained,  they  can  be  used  as  truth 
data  to  detect  and  identify  an  unknown  object.  The  following  paragraphs  provide  the  step  by  step  pro¬ 
cedure  for  training  the  algorithm  on  a  specific  object. 
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Figure  7.10:  Training  Procedure  for  Object  Detection 
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IVaining  step  1 


The  first  step  creates  8  1-D  signals,  for  i  =  {1, 2, ,  8},  using  the  Sampled  Radon  Transform 
where  A0  =  7r/8.  The  superscript  on  q^{p)  indicates  that  these  projections  are  the  original  projections 
that  are  input  to  the  first  pass  of  the  detection  algorithm.  The  training  image  is  shown  in  Figure  7.11.  It 


Figure  7.11:  Training  image  of  fighter  aircraft 


consists  of  a  fighter  aircraft  against  a  uniform  background.  Its  8  projections  are  shown  in  Figure  7. 1 2. 


IVaining  step  2 

Each  projection  is  input  to  the  matched  wavelet  algorithm  of  Chapter  5  which  creates  a  set  of  8  matched 
wavelets,  where  the  superscript  refers  to  pass  0.  The  first  step  of  the  matching  algorithm  dilates 

the  signal  such  that  there  is  a  maximum  amount  of  energy  in  the  wavelet  passband.  Because  the  first  pass 
will  primarily  find  the  low  pass  energy,  the  wavelet  matching  algorithm  is  set  such  that  the  passband  is 
narrow,  2'k/Z  <  |c<;|  <  87r/3,  which  are  the  bandlimits  for  orthonormal  wavelets.  Because  the  matched 
wavelets  from  pass  0  will  be  used  to  determine  the  rotation  of  the  object,  it  is  necessary  that  they  be 
symmetric  in  order  to  handle  the  symmetry  relationship  given  in  (7.4).  For  this  reason,  the  phase  match¬ 
ing  algorithm  will  not  be  used  in  order  to  guarantee  symmetric  matched  wavelets.  The  resultant  matched 
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wavelets  are  shown  in  Figure  7. 1 3  along  with  the  dilated  projections.  Notice  that  although  each  matched 


Figure  7.13;  Dilated  projections  and  their  associated  matched  wavelets  -  pass  0 

wavelet  has  zero  mean,  the  matching  algorithm  found  wavelets  that  behave  like  low  pass  filters,  since 
only  the  main  lobe  of  the  wavelet  overlays  the  projection.  Furthermore,  the  matched  wavelets  for  each 
projection  are  identical,  since  each  is  trying  to  detect  the  low  frequency  content.  Therefore,  the  main 
lobe  of  each  matched  wavelet  provides  information  about  the  position  of  each  projection  and  hence,  the 
object.  A  mask  function,  mi{p),  for  each  projection  is  created  by  truncating  the  matched  wavelet  out¬ 
side  its  main  lobe  and  setting  all  non-zero  values  to  1.  This  mask  function  will  be  used  later  to  filter  out 
spurious  data  not  coincident  with  the  object. 
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lyaining  step  3 

The  next  step  is  to  take  the  Wavelet  Radon  Transform,  shown  in  Figure  7.14,  using  the  matched  wavelets 
as  the  transform  operators.  As  in  Figure  7.9,  the  CWT  of  the  projection  at  0°  is  at  the  top  and  the  scale 


Figure  7. 14:  Wavelet  Radon  Transform  of  jet  aircraft  projections  -  pass  0 
parameter,  a,  increases  downward  for  each  CWT.  The  horizontal  axis  is  the  shift  parameter,  b. 

Training  step  4 

Information  regarding  the  aspect  ratio  of  the  object  can  be  extracted  from  the  maximum  scale  vector, 
which  consists  of  the  scales  at  which  the  peak  occurs  in  each  CWT  in  the  Wavelet  Radon  Transform, 
that  is, 

^max  =  {o-lnaxl^qiKnax^^)  E  {1,2, ...  ,S},  V(a,6)€5R}  (7.8) 
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where  W^{a,  h)  is  the  Wavelet  Transform  of  projection  qi{p).  Because  of  the  nature  of  the  Wavelet 
Radon  Transform,  st.max  is  scale  and  shift  invariant  and  rotation  invariant  to  8  levels.  The  maximum 
scale  vector,  amox>  for  the  fighter  aircraft  is 


28 

32 

39 


^max  — 


42 

43 


41 


37 

31 


(7.9) 


Notice  that  the  scale  moves  from  low  to  high  and  back  to  low  again.  That  is  because  the  aircraft  is  longer 
than  it  is  wide.  Projection  0  (0°)  represents  the  width  of  the  aircraft  and  projection  5  (90°)  represents 
the  length.  Because  the  width  is  smaller  than  the  length,  the  match  in  the  CWT  represented  by  the  peak 
will  occur  at  lower  a. 

At  this  point,  the  benefits  of  using  matched  wavelets  over  standard,  pre-designed  wavelets  can  be 
demonstrated.  The  peak  values  from  the  Wavelet  Radon  Transform  (WRT)  found  in  Training  step  3 
were  compared  to  those  obtained  from  the  WRT  of  the  same  projections,  but  using  Meyer’s  wavelet. 
Figure  7.15  shows  that  the  matched  wavelets  consistently  produce  higher  peaks  in  the  CWTs  of  each 
projection. 


IVaining  step  5 

Now  that  the  first  level  of  training  is  complete,  a  set  of  residual  detail  signals  can  be  constructed  from 
the  results  of  pass  0.  The  wavelet  operator  associated  with  the  peak  of  the  CWT  in  pass  0  is  used  as  truth 
for  the  information  detected  in  the  first  pass,  Assume  the  peak  in  kF*(a,  b)  occurs  at  {am,  bm)-  Then  the 
estimate  of  the  projection  information  matched  in  pass  0  is  given  as 

q^ip)  =  Wi  {am,bm)  amhi  m,(p)  (7.10) 

V.  ) 
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Figure  7.15:  Peak  values  from  the  WRT  -  matched  wavelets  vs  Meyer’s  wavelet 
where  mj(p)  is  the  mask  function  calculated  in  Step  2  above.  The  residual  detail  signals  are  given  as 

Qi{p)  =  <lUp)-Qiip)-  (7.11) 

Figure  7.16  shows  each  projection,  qi{p),  and  its  estimate,  qj  {p). 

IVaining  step  6 

The  new  set  of  residual  projection  details,  qfip)  are  input  to  the  wavelet  matching  algorithm  from  which 
a  second  set  of  matched  wavelets  is  derived,  'ipj  (p) .  In  this  pass,  however,  the  bandlimits  in  the  matching 
algorithm  are  set  to  0  <  |a;|  <  47r  and  the  phase  of  the  wavelet  is  set  to  the  phase  of  its  corresponding 
projection.  The  matching  algorithm  will  therefore,  generate  asymmetric,  dyadic  wavelets,  which  is  ap¬ 
propriate  since  the  residual  projection  detail  signals  should  now  contain  predominantly  high  frequency 
content.  Figure  7. 17  shows  the  matched  wavelets  (solid  line)  and  their  corresponding  projections,  dilated 
for  maximum  energy  (dashed  line). 
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Figure  7.16:  Projections  and  their  estimate  -  pass  0 
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Figure  7.17:  Dilated  projections  and  their  associated  matched  wavelets  -  pass  I 
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IVaining  step  7-8 


Steps  3  and  5  are  repeated  to  obtain  g?(p)  and  qf{p),  which  are  the  projection  estimate  of  pass  2  (Fig¬ 
ure  7.18)  and  the  new  projections  for  pass  3  (if  needed).  Figure  7.18  shows  the  projection  detail  signals 
from  pass  1  (solid  line)  and  the  estimates  of  the  signal  matched  (dashed  line).  The  data  generated  in  this 


Figure  7. 1 8;  Projections  and  their  estimate  -  pass  1 

section  can  now  be  applied  to  an  unknown  object  in  order  to  determine  if  that  object  is  the  fighter  aircraft 
used  for  training. 

7.4.2  Object  Detection 

The  algorithm  for  detecting  and  identifying  an  object  developed  in  this  section  is  a  modification  of  the 
backprojection  image  reconstruction  algorithm  described  in  Section  7.2.2,  and  the  Wavelet  Radon  Trans¬ 
form  developed  in  Section  7.3.  Figure  7.19  shows  the  algorithm  developed  in  this  section.  Notice  its  sim- 
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Figure  7. 1 9:  Object  Detection  Algorithm 
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ilarity  to  the  backprojection  image  reconstruction  algorithm  shown  in  Figure  7.5.  The  key  difference  is 
that  the  |7y|  filter  is  replaced  by  matched  wavelet  detection.  The  following  paragraphs  will  explain  each 
of  the  steps  in  the  object  detection  algorithm  using  as  an  example  image,  the  same  image  ( Figure  7.11) 
used  for  training. 

Detection  step  1 

The  first  step  constructs  8  projections  from  the  object  image.  Since  the  test  image  is  the  same  as  that 
used  in  training,  the  projections  will  obviously  be  the  same  (Figure  7.12). 

Detection  step  2 

Step  2  calculates  the  CWT,  Wq{a,  b),  of  the  projections  using  the  matched  wavelets  for  pass  0  found 
during  training.  Again,  since  the  example  image  in  this  case  is  identical  to  the  training  image,  the  CWT 
will  be  that  of  Figure  7. 14. 

Detection  step  3 

The  peak  scale  vector,  defined  in  (7.9),  is  derived  from  the  CWT  and  will  be  used  to  determine  whether 
the  object  is  the  same  as  the  training  object,  and  if  so,  its  rotation  and  scale.  Let  a  be  the  peak  scale 
vector  corresponding  to  the  CWT  found  in  this  step.  If  the  test  object  is  a  rotated  and  scaled  version  of 
the  training  object,  then  a  will  be  a  rotated  and  scaled  version  of  a.  So,  in  order  to  determine  rotation, 
a  is  circularly  correlated  with  a  found  in  Training  Step  4,  and  then  normalized  by  the  maximum  of  the 
result,  that  is, 

a*  a 

s=  - p - r  (7.12) 

sup  [a  *  aj 

where  s  =  {s(()|0  <  s(i)  <1,  i  =  {1, 2, ...  ,8}}  The  angle  of  rotation  will  correspond  to  the 
location  of  the  1  (maximum  value)  in  the  circular  correlation.  For  the  example  test  object,  a  =  a  = 
{28, 32, 39, 42, 43, 41, 37, 31}.  The  normalized  circular  correlation,  shown  in  Figure  7.20,  indicates  that 
the  proposed  angle  of  rotation  is  0°,  as  it  should  be.  The  angle  is  considered  a  proposed  angle  because 
it  has  not  yet  been  determined  that  the  object  in  question  is  in  fact  the  training  object. 
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Figure  7.20:  Normalized  Circular  Correlation,  s  -  0° 


Once  the  proposed  angle  of  rotation  is  determined,  a  can  be  circularly  reordered  to  align  with  a.  In 
other  words,  if  the  proposed  rotation  were  22.5°,  then  0.2  would  correspond  to  ai,  03  to  02  and  so  on, 
finishing  with  oi  corresponding  to  as,  due  to  the  symmetry  of  the  Radon  Transform.  The  reordered  a, 
call  it  a,  can  be  used  to  find  the  proposed  scale  of  the  test  object  and  a  decision  criteria  for  detection.  If 
the  test  object  is  a  scaled  and  rotated  version  of  the  training  object,  a,  having  been  corrected  for  rotation, 
will  be  a  scaled  version  of  a  and  the  scale  factor  across  i  will  be  a  constant  in  the  ideal  case.  Let  7  be 
defined  as 

7  =  {7i|7i  =  ai/oi,  i  =  {1,2,...,8}}  (7.13) 

and  let 

1  ® 

(7.14) 
and 


8^, 

Z— 1 


(7.15) 

be  the  mean  and  standard  deviation  of  7.  If  a  is  below  some  threshold,  then  the  test  object  is  detected, 
that  is,  it  is  the  same  as  the  training  object,  and  ^  is  the  scale  factor  of  the  test  object  to  the  training  object. 
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Detection  step  4 


Once  the  test  object  is  determined  to  be  the  training  object,  the  remainder  of  the  algorithm  can  provide 
additional  information.  Assume  the  peak  of  W*(o,  b)  occurs  at  {a!^,  b\^).  The  estimate  of  the  projec¬ 
tion  information,  q^{p),  calculated  using  (7.10),  where  mi{p)  is  the  mask  function  described  in  Training 
step  5,  is  the  detected  information.  The  residual  details  signals,  ql{p),  given  by  (7.11)  represent  the  de¬ 
tails  remaining  in  the  image  after  pass  0.  Each  of  these  sets  of  projections  are  backprojected  creating 
a  detected  image,  /<i(x,y),  and  a  residual  detail  image,  f^{x,y).  Figure  7.21a  shows  the  detected  im¬ 
age  /d(x,  y)  overlayed  on  the  test  image.  Notice  that  the  detected  image  reflects  the  position,  aspect 
ratio,  rotation,  and  scale  of  the  target.  Figure  7.21b  shows  the  resultant  detail  image  after  the  first  pass. 
Both  backprojected  images  were  thresholded  at  70%  of  their  maximum  to  eliminate  the  noise  inherent 
to  backprojection. 

Detection  step  5 

Before  repeating  steps  3-4  for  qjip),  the  matched  wavelets,  V’i(p)v  calculated  in  Training  step  6  have 
to  be  adjusted  for  rotation.  First,  they  have  to  be  reordered  in  the  same  way  that  a  was  reordered  so 
that  they  will  align  with  the  test  projections.  Second,  some  of  the  matched  wavelets  must  be  flipped 
along  their  horizontal  axis.  The  training  projections  were  collected  at  0  <  P  <  tt.  If  the  test  object  is 
rotated  an  angle,  (p,  relative  to  the  training  image,  then  it  is  as  though  the  test  projections  were  collected 
at  (f  <  6  <  +  (f.  Because  g^+;r(p)  =  g^(-p),  the  test  projections  collected  for  0  >  tt  are /9- 

reversed.  Therefore,  their  corresponding  matched  wavelet  must  be  p-reversed  so  that  they  will  match 
their  corresponding  projections  in  the  CWT.  This  procedure  was  not  necessary  in  the  first  pass  because 
the  matched  wavelets  were  symmetric. 

Once  the  matched  wavelets  for  the  second  pass  are  adjusted.  Detection  steps  3-4  are  repeated  with 
q}{p)  as  input.  The  backprojection  results  are  shown  in  Figure  7.21c  and  d.  The  results  in  Figure  7.21 
represent  a  template  for  a  human  observer  to  use  to  further  identify  a  test  object. 
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7.4.3  Detection  Results  for  an  Unknown  Object 

To  test  all  the  functions  developed  in  the  previous  section,  a  test  image,  shown  in  Figure  7.22,  was  con¬ 
structed  by  rotating,  scaling,  shifting  and  adding  noise  to  the  fighter  aircraft  shown  in  Figure  7.11.  The 


Figure  7.22:  Test  image  for  detection 

angle  of  rotation  is  1 12.5° ,  and  the  scale  factor  is  1 .5.  The  aircraft  was  translated  15  pixels  up  and  30  pix¬ 
els  to  the  right  of  center  and  23dB  of  uniform  white  noise  was  added  making  the  calculated  SNR=  2.1dB. 
Figure  7.23  shows  the  test  projections  and  their  detected  estimate  found  from  Detection  steps  2-4.  The 
Wavelet  Radon  Transform  from  pass  0  is  shown  in  Figure  7.24.  The  same  characteristic  sinusoidal  path 
seen  in  the  Radon  Transform  of  a  shifted  point  is  evident  in  the  WRT  of  the  translated  object. 
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Figure  7.24:  Wavelet  Radon  Transform  of  Test  image  projections  -  pass  0 


The  maximum  scale  vector  for  the  test  image  is  given  as 


64 

64 

60 

56 

48 

41 

54 

60 


(7.16) 


The  results  of  circularly  correlating  a  with  a  in  (7.9)  is  shown  in  Figure  7.25.  The  peak  corresponds  to 


Figure  7.25:  Normalized  Circular  Correlation,  s  -  Test  image 
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6  =  112.5°,  as  expected.  The  vector,  7,  is  calculated  using  (7.13)  to  be 


1.464 

1.688 

1.538 

1.524 

1.488 

1.463 

1.514 

1.548 


(7.17) 


The  mean,  //  =  1.529,  and  standard  deviation,  a  =  0.0671,  of  7  give  the  scale  factor  and  goodness  of  fit 
of  the  test  object  to  the  training  object.  The  mean  is  close  to  the  actual  scale  factor  of  1.5  and  a  is  very 
small  indicating  that  the  object  is  in  fact  the  fighter  aircraft  used  for  training. 

The  results  of  backprojecting  both  the  detected  signals,  qf  {p),  and  the  detail  signals,  q\  {p)  are  shown 
in  Figure  7.26a  and  b.  Notice  that  Figure  7.26b  looks  very  much  like  Figure  7.21b.  The  estimated  details 
survived  scaling,  rotation,  translation  and  noise.  These  details  can  be  used  as  additional  templates  to  a 
human  observer  to  verify  the  automatic  detection  results  or  for  further  more  detailed  classfication. 

The  WRT  of  q]  (p)  using  the  matched  wavelets  for  pass  1  is  shown  in  Figure  7.27.  The  projection  de¬ 
tails  and  their  detected  estimates  are  shown  in  Figure  7.28.  The  backprojected  results  of  detection  in  pass 
1  and  the  new  details  are  shown  in  Figure  7.26  c  and  d.  The  results  of  pass  1  are  shown  in  Figures  7.28, 
7.27,  and  7.26c  and  d. 
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Figure  7.26:  Results  of  backprojected  detection  and  detail  signals  -  Test  image 


Figure  7.27:  Wavelet  Radon  Transform  of  Test  image  projections  -  pass  1 


137 


Figure  7.28:  Test  image  projections  and  their  estimate  -  pass  1 
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Chapter  8 


Summary 


The  research  results  presented  in  this  dissertation  show  by  proof  and  example  that  it  is  possible  to  design 
bandlimited  wavelets  that  optimally  match  a  signal  of  interest  in  the  least  squares  sense.  The  algorithm 
developed  can  be  constrained  such  that  the  resultant  wavelets  either  form  an  orthonormal  basis  of  (3?) 
or  are  dyadic  wavelets  that  satisfy  the  stability  condition.  The  development  of  the  wavelet  matching 
algorithm  provides  insight  into  the  characteristics  of  both  the  amplitude  and  phase  of  an  orthonormal 
wavelet  spectrum.  The  wavelet  matching  algorithm  is  superior  to  existing  wavelet  design  techniques 
in  that  it  places  the  design  emphasis  on  the  shape  of  the  wavelet  and  the  RMS  error  with  respect  to  the 
signal  of  interest.  The  target  detection  and  identification  algorithm  was  developed  to  be  shift,  scale  and 
rotation  invariant  with  optimal  detection  capability  by  combining  the  Radon  Transform,  backprojection 
and  the  Continuous  Wavelet  Transform  with  matched  wavelets.  It  was  demonstrated  to  work  well  with 
a  target  against  a  benign  background  embedded  in  white  noise.  The  residual  projection  details  were 
shown  to  survive  rotation,  shift,  scale  and  noise.  The  second  level  of  residual  projections,  however,  did 
not  survive,  which  was  due  to  their  small  amplitudes  relative  to  the  projections  from  the  first  pass.  The 
matched  wavelets  increased  the  detection  performance  as  expected  since  they  were  specifically  designed 
to  match  the  projections  under  analysis.  As  the  targets  and  backgrounds  become  more  complicated,  this 
optimal  performance  should  prove  beneficial  in  detecting  targets  in  clutter. 
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8.1  Future  Research 


The  next  logical  development  step  for  the  wavelet  matching  algorithm  is  to  extend  it  to  bi-orthogonal 
wavelets,  that  is,  wavelets  that  reside  in  a  multiresolution  analysis  (MRA)  context  whose  duals  reside 
in  a  dual  MRA.  The  QMF  constraints  exploited  in  the  development  of  the  matching  algorithm  in  this 
dissertation  become  much  more  complicated  for  bi-orthogonal  wavelets.  This  added  complication  will 
need  to  be  overcome  before  the  algorithm  can  be  extended  to  bi-orthogonal  wavelets. 

The  target  detection  and  identification  algorithm  can  be  improved  to  include  targets  in  clutter.  Sur¬ 
rounding  clutter  will  add  significant  unwanted  energy  to  the  projections.  It  may  be  necessary  to  take 
windowed  Radon  Transforms  to  eliminate  the  unwanted  energy  from  other  parts  of  the  scene.  The  iden¬ 
tification  capability  based  on  aspect  ratio  can  be  improved  to  full  classification  capability  by  applying 
pattern  classification  techniques  to  the  residual  detail  signals  from  the  first  pass.  Also,  increasing  the 
number  of  projections,  while  increasing  the  computational  burden,  may  increase  the  algorithm’s  perfor¬ 
mance  significantly. 

8.2  Conclusion 

Wavelet  theory  and  its  applications  are  exciting  and  growing  fields.  The  zoom  in,  zoom  out  capability 
of  the  wavelet  transform  and  the  ability  to  build  matched  wavelets  makes  their  application  to  image  pro¬ 
cessing  even  more  significant.  The  engineering  community  is  still  on  the  threshold  of  developing  solid, 
useful  applications  of  wavelets  in  signal  processing  that  truly  take  advantage  of  all  that  wavelets  have  to 
offer.  Like  many  other  engineers  and  scientists,  much  of  the  research  conducted  in  support  of  this  dis¬ 
sertation  was  focused  on  developing  an  in-depth  understanding  of  wavelets  and  their  characteristics.  It 
seems  clear  that  once  the  engineering  community  has  the  same  depth  of  understanding  of  wavelet  theory 
as  they  do  for  Fourier  theory,  they  will  develop  new  applications  that  are  ideally  suited  to  wavelets  and 
the  advantages  they  hold. 
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Appendices 


Appendix  A:  Proof  -  Guaranteeing  Orthonormality 


Substituting  lj  =  tj  +  lirm  in  (5.37)  gives 

1  for  7/  +  27rm  =  0 

(_  I'®’  +  27rm))  for  t]  +  2nm  7^  0 

Summing  both  sides  of  (A-1 )  over  m  gives 


|$(7?  +  27rm)p  =  ^ 


(A-1) 


^  \^{ri  +  2'Km)\^  =\ 


(A-2) 


1  r}  +  2‘nm  —  0 

i  EmEn  1’^  (2""''H»7  +  27rm))|^  r/-|-27rm/0 
The  left  side  of  (A-2)  is  the  Poisson  summation,  which  by  (5.16)  must  is  1  everywhere  when  4>j,k  is  an 
orthonormal  basis  of  Vj.  Therefore,  because  of  the  equality  in  (5.37),  (f)j^k  is  an  orthonormal  basis  of  Vj 
if  and  only  if 

00  00  2 

^  ^|^r(2"+i(^  +  27rm))|  =1.  (A-3) 

n 

□ 


m=~oo  n—0 


Appendix  B:  Proof-Bandlimited  $ 

The  Fourier  transform  of  the  scaling  function’s  2-scale  relation  is  given  as 

$M  =  h(|)*(|).  (B.1) 
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From  this  expression,  it  is  clear  that  the  maximum  bandwidth  of  $(a;)  is  equal  to  the  maximum  band¬ 
width  of  one  period  of  H{u)/2).  Let  the  maximum  half-bandwidth  of  be  given  as  q^tt  then  the 
maximum  half-bandwidth  of  one  period  of  H{u))  is  also  aix.  The  maximum  half-bandwidth  of  $(a;/2) 
is  2a7r.  In  order  for  (B-1)  to  be  satisfied  for  bandlimited  scaling  functions,  ^{uj/2)  must  not  extend  into 
the  second  period  of  Ff  (a;/2),  which  implies 


2a7r  <  in  —  an 

“  -  i 

The  maximum  scaling  function  bandwidth  is  therefore  |a;|  <  47r/3. 

□ 


(B-2) 


Appendix  C:  Proof-Bandlimited  ^ 


‘^max  and  be  the  upper  frequency  bounds  of  bandlimited  <4^  and  respectively.  From  (5.37) 
|$(u;)|2  =  |^(2a;)|2  +  |5r(4a,)|2  +  ...  for  a;  7^  0.  (C-1) 


Clearly,  the  upper  frequency  bound  of  $(a;)  is  due  to  ^{2uj)  only.  Therefore, 


ijr  _  „  $  _  StT 

^max  ^^max  n  • 


If  (p{x  —  k)  is  orthonormal,  its  Poisson  summation  must  be  1  everywhere. 


00 

|$(w  +  27rm)p  =  1 

m=-oo 

which  for  bandlimited  scaling  functions  (Theorem  6)  implies 


(C-2) 


(C-3) 


|$(a;)|  =  l  for^<|u;|<f.  (C-4) 

Let  the  lower  frequency  bound  of  ’4'  be  given  as  so  that  the  bandlimits  of  ^>(0;)  is  <  |a;|  < 
87r/3.  Let  $1(0;)  be  the  portion  of  $  between  7r/3  and  27r/3.  From  (C-1)  and  (C-2) 

|$i(w)|2  =  \^{2u)\‘^  +  |«^(4u;)|2  =  1  for  I  <  |u;|  <  f .  (C-5) 
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Let  $2(^4')  be  the  portion  of  $  between  tt/C  and  7r/3  and  the  portion  of  'S'  between  tt/S  and  27r/3, 
then 

|$2HP  =  \^{2u)\^  +  |'S'(4a;)|2  +  \^{8u})\‘^  =  1 

=  |'S'2(2a;)P  +  |$i(2a;)|2  ^1  forf<|a;|<f.  (C-6) 

=  |'S'2(2u;)P +  1  =1 

So, 


|'®'2MP  =  0  forf<|u;|<f.  (C-7) 

We  can  generalize  this  approach  by  letting  $j(w)  be  the  portion  of  $  between  2tt/3  ■  2^  and  47r/3  •  2^ 
and  'S'j(w)  be  the  portion  of  'S'  between  47r/3  •  2^  and  87r/3  ■  2^,  for  j  =  {2, 3, . . . ,  00},  then 


=  Ei=o  l»(2‘+'u.)p  =  1  for  ^  <  li^l  <  ^. 
However,  for  any  value  of  j,  has  already  been  shown  to  be  0,  so  that 

A;=0  n=2 

=  |^j(2u;)|2  +  1. 


(C-8) 


(C-9) 


Substituting  back  into  (C-8)  gives 


=  |«',(2u;)p  -M  =  1  for  <  |u;|  <  ^.  (C-10) 

Therefore, 

^j{2ij)  =0  for  ^  <  |u;|  <  ^  (C-1 1) 

and 

=  0  for  ^  <  |tj|  <  (C-12) 

Combining  all  for  {y  =  2, 3, . . . ,  00}  gives 

|’®'(a;)|2  =  0  for0<|w|<f.  (C-13) 

Therefore,  =  27r/3  and  the  bandlimits  on  'I'  for  bandlimited  $  in  an  orthonormal  MRA  is 

27r  I  ^  Stt 

y  <  kl  <  y-  (C-14) 

□ 
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Appendix  D:  Proof-Guaranteeing  an  Orthonormal  MRA 


Let  |^(A:Aa;)|  be  a  sampled  wavelet  spectram  in  an  orthonormal  MRA  with  sample  spacing  Aa;^  = 
27r/2^.  According  to  Theorem  4,  we  can  calculate  the  scaling  function  spectrum  at  any  value, 
|$(fc7r/2^)|,  using  (5.29),  where  £  =  {0, 1, ,  M}.  Let  Aa;  =  7r/2^  in  (5.29),  then 

f  2^‘''^raAw\  P 


|$(nAw)p  =  ^ 

p=0 


\  2^  J 

Substituting  nAw  =  nAw  +  27rm  and  summing  over  m  gives 

'2^+1 


(D-1) 


|$(nAa;  +  27rm)p  =  Y  S 

ri=  — OO  m=  — oop=0 


2P 


{nAu  +  27rm) 


(D-2) 


The  left  side  of  (D-2)  is  the  Poisson  summation  sampled  at  Aa;  and  is  therefore  equal  to  1 .  Substituting 
back  in  for  Aw  gives 


m=— oop=0  ^  ^ 

Let  y(fc)  =  |'^(fcAw$)P  where  Aw$  =  2tt/2^  ,  then  the  constraint  on  Y  becomes 

oo  e  /^M  . 

E  Ei^  |r("  +  2'«m)  =1 

Tn~—oopz=0  \  J 

for^  =  {0,l,...,M}. 

Rewriting  (5.42)  in  Corollary  7  for  sampled  wavelet  spectra  gives 

27r  .  I  ^  Stt 

T  -  ^  T- 

Dividing  through  by  Aw  =  27r/2^  gives  the  index  constraint  on  k. 

oM  9M+2 

Substituting  the  index  of  Y  in  (D-4)  gives  the  bandlimit  constraint  on  the  samples  of  Y. 


=  1. 


2M 


< 


2M 


(n  + 


< 


2M+2 


(D-3) 


(D-4) 

(D-5) 


(D-6) 


(D-7) 


(D-8) 


□ 
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Appendix  E:  Proof-Matched  Wavelet  Amplitude 


The  problem  of  finding  the  optimal  Y  given  the  error 

(W  -  Y)^(W  -  Y) 


E  = 


and  the  equality  constraints 


W^W 


AY  =  1 


can  be  solved  easily  using  Lagrangian  multipliers.  The  Lagrangian  is  given  as 


L  = 


(W  -  Yf(W  -  Y) 


+  A^(AY  -  1). 


W^W 

The  error  function,  E,  is  minimized  by  setting  VyL  =  0,  VxL  =  0  and  solving  for  Y  [21]. 

2 


VyL  = 


(W  -  Y)  +  A^  A  =  0. 


W^W 

Setting  VxL  =  0,  we  get  back  our  equality  constraints. 

AY  =  1. 

Multiplying  (E-4)  by  A  and  substituting  (E-5)  gives  a  solution  for  A. 

^  (AW-AY)  +  AA^A  -  0 


W^W 

2 


W^w 

2 


(AW-1)  +  AA^A  =  0 


(AA^)-^(l  -  AW)  =  A. 


Substituting  (E-6)  back  into  (E-4)  gives  the  solution  for  Y. 


(E-1) 


(E-2) 


(E-3) 


(E-4) 


(E-5) 


(E-6) 


Y  =  W  -I-  A'^  (AA^)-i(l  -  AW) 


(E-7) 


and  substituting  (E-7)  into  (E-1)  gives  the  error  in  the  match, 

^  _  -(A^(AA^)-i(l- AW))^(-l)(A^(AA^)-i(l- AW)) 

W^w 

_  (1  -  AW)^(AA^)-^AA^(AA^)-M1  -  AW) 

W^W 

_  (1  -  AW)^(AA^)-1(1  -  AW) 

W^W 


(E-8) 
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Because  Y  is  a  power  spectrum,  we  must  add  the  inequality  constraints 


Y{k)  >0  for  all  k. 


(E-9) 


Note:  If  an  element  ofY  is  negative,  then  we  must  add  an  equality  constraint  forcing  that  element  to  0. 
This  is  done  by  appending  a  0  to  1  and  adding  a  row  to  the  bottom  of  A  containing  all  O’s  and  a  1  in  the 
column  corresponding  to  the  negative  element. 

Notice  in  (E-7)  and  (E-8)  that  Y  and  E  are  very  much  dependent  on  the  amplitude  of  W.  Since  the 
matched  wavelet  must  match  the  desired  signal  up  to  a  scale  factor,  K,  as  shown  in  (5.4),  we  can  scale 
W  so  that  E  is  minimized.  Let  the  scale  factor  be  1/a  so  that  the  desired  signal  spectrum  is  W/a.  The 
error  function  becomes 


E{a)  = 


(1  -  -  iAW) 


and  the  solution  for  Y  becomes 


(E-10) 


Y  =  -  W  +  A’^(AA^)-1(1  -  -AW). 
a  a 


(E-11) 


Setting  dE{a)/da  =  0  gives 


d  (1  -  iAW)^(AA^)-i(l  -  iAW) 
da  4W^W 

d  (al  -  AW)^(AA^)-Hal  -  AW) 


da 


W^W 


=  0 


=  0 


[(al-AW)^(AA^)-4  +  l^(AA^)-i(al-AW)]  =  0 

[l^(AA^)~^(al- AW)]  =  0.  (E-12) 


W^W 

Solving  for  a  in  (E-12)  gives  the  expression  for  the  scale  factor  that  will  minimize  the  error  function  E, 

l^(AA’^)-iAW 


a  = 


(E-13) 


lT{AAT)-il  • 

Because,  the  conditions  on  from  Theorem  8  are  satisfied  through  the  constraints  of  the  optimization 
equations,  the  scaling  function  spectrum  calculated  by  way  of  (5.29)  generates  an  orthonormal  MRA. 
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Appendix  F:  Proof-Properties  of  A$,  A$,  and  A 

The  proofs  for  (5.62),  (5.63)  and  (5.64)  are  trivial.  Since  ^(a;),  ipix),  and  hk  are  all  real  quantities,  their 
Fourier  Transforms  are  Hermitian[3],  that  is,  their  magnitude  is  even  and  their  phase  is  odd.  Therefore, 
and  9h  are  odd  functions.  Since  the  derivative  of  an  odd  function  is  an  even  function  if  follows 
that  A$(ci;),  A^(u;),  and  A(a;)  are  all  even  functions.  In  order  to  prove  (5.65)  and  (5.66),  let  Ao(w)  = 


A(u;)  —  A,  where 


so  that 


1  r 


/X{oj)du 

-TT 


1 


The  average  value  of  A$(a;)  is  given  by 


/OO  fOO  ^ 

A^{ij)duj  =  / 

roo  OO 

=  /  S 


dw  =  0. 


4“  A  duj 


=  P-MM 


A  ^2- 


=  A. 


Similarly  for  A^  (w) 

/OO 

A\if{u;)du; 

-OO 


27r7_^  2  2  V2  /  ^2 

=  (^“  (I  +  +  P""  (•'"  (I:) 


11 
2  2 


A  +  A  ^  2" 

771=2 


1  1^  1^ 
~2  2"^  ^  2^ 
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Appendix  G:  Continuous  Wavelet  Transform  Implementation 

The  continuous  wavelet  transform  is  given  by  the  following  expression, 

Wf{a,b)=j  f{x)a~^'tp  (G-1) 

The  wavelets  used  in  this  dissertation  are  discrete  functions,  'ip{k)  =  ip{kAx),  calculated  using  the 
matching  algorithm  of  Chapter  5.  Implementing  the  CWT  with  discrete  quantities  is  done  with  sum¬ 
mation  rather  than  integration  and  the  result  is  actually  only  an  approximation  to  the  CWT.  The  approx¬ 
imation  to  (G-1)  is  given  by 

Wf{m,n)=  f{kAx)'ip\^ — — j  (G-2) 

where  N  is  the  number  of  points  in  f{kAx)  and  'ij}{kAx).  In  order  to  simplify  (G-2),  assume  Ab  =  Ax, 
which  provides  the  finest  translation  steps  possible,  since  the  shift  parameter  cannot  step  any  smaller  than 
Ax.  Rewriting  (G-2)  and  dropping  the  Ax’s  gives 

N/2-1 

k=-N/2  ^  ' 

which  can  be  written  in  matrix  notation  as 

Wfim)  =  (G-4) 

where  Wf  (m)  is  row  m  of  the  approximate  CWT  and  has  length,  iV,  f  is  a  iV  x  1  vector  consisting  of  the 
samples  of  the  input  signal,  f{x),  and  isaN  x  N  matrix  where  each  row  consists  of  the  sampled 
wavelet,  dilated  to  scale  mAa  and  shifted  by  n,  that  is,  -  n) /mAa. 

All  approximate  CWTs  in  this  dissertation  will  be  calculated  for  64  discrete  scales.  Furthermore,  as 
in  Chapters  6  and  7,  JV  =  512,  Arc  =  1/32,  and  Aw  =  tt/S.  Let  V’(fc)  be  the  discrete  wavelets  calculated 
by  the  matching  algorithm  of  Chapter  5  with  the  parameters  given  above.  Because  Arc  =  1  /32,  ^(/s)  can 
be  decimated  up  to  5  times  giving  \/^ip{32k),  where  now  Ax  —  1.  In  this  implementation,  however, 
'ip{k)  is  decimated  only  4  times  since  a  fifth  decimation  produces  a  function  with  only  one  non-zero  value. 
So,  if  ^/l6'^p{16k)  is  used  to  calculate  the  first  row  of  the  CWT  which  would  represent  the  smallest  scale, 
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then  \/8^(8A;)  will  be  used  to  calculate  the  second  row,  \/4^(4A;)  the  fourth  row,  y/2'ip[2k)  the  eighth 
row  and  the  sixteenth  row.  Upsampling  and  interpolating  gives  'il){k/2)^  which  will  be  used 
to  calculate  the  32nd  row  of  the  CWT  and  upsampling  and  interpolating  again  gives  '^(A:/4),  which  will 
be  used  to  calculate  row  64  of  the  CWT.  There  are  now  discrete  wavelets  that  will  calculate  the  CWT  at 
rows  2^  for  j  =  0, 1, . . . ,  6.  The  wavelets  that  will  be  used  to  calculate  the  other,  non-dyadic  rows,  are 
found  by  interpolating  the  wavelets  at  the  dyadic  rows.  For  instance,  to  find  the  wavelet  corresponding 
to  the  23rd  row  of  the  CWT,  the  wavelet  corresponding  to  row  16  is  upsampled  by  a  factor  of  23/16, 
which  is  done  by  upsampling  and  interpolating  by  a  factor  of  23/1  and  taking  every  16th  sample.  Now 
that  the  wavelets  corresponding  to  every  scale  in  the  output  CWT  has  been  calculated,  they  can  be  used 
in  (G-4)  to  find  the  approximate  CWT. 
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